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1 . INTRODUCTION 

1. 1 Objective  and  Scope 

The  main  purpose  of  this  investigation  is 
to  develop  a numerical  procedure  for  accurately  solving 
two-point  boundary  value  problems  using  an  initial  value 
or  "shooting"  approach  (Roberts,  1962,  and  Salvadori  and 
Baron,  1961)*  when  rapidly  growing  solutions  are  present. 

The  method  is  applicable  to  a wide  variety  of  practical 
engineering  problems,  but  herein  will  be  demonstrated  only 
in  the  solution  of  a fourth-order  system  of  differential 
equations  with  two  growing  solutions  in  each  direction,  and 
in  the  determination  of  natural  frequencies  and  mode  shapes 
of  sperical  shell  segments. 

It  is  well  known  that  in  the  application  of  a 
straightforward  shooting  method,  difficulties  are  encoun- 
tered when  rapidly  growing  solutions  are  present.  Such 
rapidly  growing  solutions  always  exist  when  there  is  an 
edge  effect,  such  as  in  long,  slender  beams  on  elastic 
foundations  or  in  the  flexural  theory  of  thin  shells. 

Several  methods  have  been  proposed  to  treat  cases  involving 
rapidly  growing  solutions,  none  of  them  entirely  satisfactory. 


Names  and  dates  in  parenthesis  refer  to  entries  in  the 
List  of  References 
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What  is  proposed  here  is  a new  method  for 
initial  value  solution  of  ordinary  differential  equations 
which  have  edge  effects  in  the  solution.  In  essence,  a 
rapidly  growing  solution  of  a differential  equation  or 
system  of  differential  equations  is  first  found.  Then  the 
differential  equation  is  altered  with  the  aid  of  the  rapidly 
growing  solution,  so  that  the  remaining  differential  equation 
does  not  contain  this  solution.  Since  the  difficulty  arising 
with  a rapidly  growing  solution  lies  in  the  fact  that  the 
growing  solution  dominates  all  shooting  solutions,  a method 
which  prevents  that  from  occurring  is  an  attractive  way  to 
eliminate  the  difficulty. 


1.2  Initial  Value  Solutions  of  Two-Point  Boundary  Value 
Problems 


Two-point  boundary  value  problems  occurring  in 
various  areas  of  engineering,  physics  and  applied  mathe- 
matics generally  cannot  be  solved  analytically.  Therefore, 
it  becomes  necessary  to  use  numerical  methods  to  find  the 
solution.  Some  of  the  numerical  methods  being  used  include 
the  following:  variational  methods,  the  method  of  colloca- 
tion, and  the  finite  difference  method  (Roberts,  1972). 

The  application  of  these  methods,  however,  normally  require 
large  amounts  of  computer  time  and  storage  space. 
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An  alternative  way  to  solve  two-point  boundary 
value  problems  is  to  replace  them  with  a set  of  initial 
value  problems.  The  procedure  merely  involves  numerically 
integrating  the  differential  equation  from  the  initial 
boundary  to  the  final  boundary  a number  of  times.  The 
integrations  start  with  given  and  assumed  independent 
initial  data/  the  total  number  of  initial  conditions  always 
being  equal  to  the  sum  of  the  number  of  given  initial  and 
the  number  of  given  final  boundary  conditions.  The  solution 
to  the  boundary  value  problem  in  question  will  then  be  ob- 
tained either  by  successive  approximations  to  the  unknown 
initial  data  or  by  combining  the  solutions  solving  simul- 
taneous equations.  This  method  is  referred  to  as  a shooting 
method  (Roberts,  1972).  Solutions  obtained  by  this  method 
do  not  present  difficulties  as  long  as  no  rapidly  growing 
solutions  exist.  Whei;  rapidly  growing  solutions  are  present, 
the  application  of  the  shooting  method  as  described  above  is 
ineffective.  The  solutions  lose  their  independence  because 
the  rapidly  growing  part  of  the  solution  eventually  dominates 
in  every  initial  value  solution.  Consequently,  since  the 
solutions  are  not  effectively  independent,  from  a numerical 
point  of  view,  they  cannot  be  accurately  combined  to  form 
the  correct  solution. 

Variants  of  the  shooting  method  have  been  devised 
to  avoid  the  problem  caused  by  rapidly  growing  solutions. 
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Among  these  are  the  subinterva1  technique,  the  suppression 
technique,  the  field  method,  and  the  invariant  imbedding 
method. 

Cohen  (Cohen,  1964)  and  Kalnins  (Kalnins,  1964) 
developed  methods  whereby  the  total  interval  is  divided 
into  a number  of  subintervals,  over  each  of  which  initial 
value  problems  are  solved  separately  using  the  same  system 
of  differential  equations.  After  the  integrations  have 
been  performed  over  each  subinterval,  continuity  conditions 
are  imposed  on  appropriate  quantities  at  the  ends  of  each 
segment,  thereby  forming  a simultaneous  system  of  linear 
equations . 

The  suppression  method  developed  by  Zarghamee  and 
Robinson  (Zarghamee,  1967)*  re-orthogonalizes  the  solutions 
whenever  they  grow  beyond  a certain  predetermined  limit. 

As  in  the  ordinary  shooting  method,  solutions  are  obtained 
corresponding  to  each  set  of  initial  conditions  by  numerical 
integration.  For  example,  let 

IV  4 

w - 0 W = 0 (1.1) 

be  the  governing  differential  equation.  Further  let 

* 

The  suppression  method  was  also  developed  independently  by 
Goldberg  (Goldberg,  1965). 
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W (0)  = WQ 

W' (0)  = 0 

be  the  actual  initial  conditions  and 

W' ' (0)  =0 

W"  ' (0)  = 1 

and  W’ ' (0)  =1 

W‘  " (0)  = 0 

be  two  assumed  sets  of  independent  initial  conditions.  Two 
solutions  are  obtained,  one  for  each  set  of  assumed  conditions, 
but  both  satisfying  the  actual  initial  conditions.  Integration 
of  equation  (1.1)  is  started  at  the  initial  boundary  and  contin- 
ued until  the  solutions  grow  beyond  some  predetermined  limit. 

When  this  happens,  a point  of  suppression  is  introduced.  Two 
new  sets  of  conditions  are  assumed.  The  two  independent  solu- 
tions are  then  linearly  combined  to  satisfy  each  set  of  new 
conditions,  thereby,  forming  two  new  independent  solutions. 

These  solutions  are  propagated  until  another  point  of  suppres- 
sion is  introduced  and  the  solutions  are  again  re-orthogonalized. 
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After  several  suppressions,  the  remaining  solutions  will 
be  very  close  to  the  correct  solution  which  does  not  contain 
the  growing  function.  It  is  obvious  that  in  this  method 
fewer  solutions  must  be  carried  along  than  in  the  Cohen  or 
Kalnins  method. 

The  field  method  (Miller,  1967,  and  Berezin  and 
Zhidkov,  1960)  uses  a technique  whereby  a differential 
equation  which  has  a growing  solution  is  replaced  by  a 
set  of  first-order  differential  equations,  involving  the 
field  functions,  which  do  not  possess  growing  solutions. 

A linear  differential  equation  of  order  2n  with  n boundary 
conditions  given  at  a and  n at  b is  first  replaced  by  a 
system  of  n second-order  differential  equations.  This  is 
accomplished  by  substituting  new  variables  for  the  dependent 
variable  and  its  derivatives.  The  system  of  n second-order 
differential  equations  is  in  turn,  replaced  by  a system  of 
n(n  + 1)  first-order  differential  equations  which  must  be 
integrated  from  a to  b and  a second  system  of  n first-order 
equations  which  have  to  be  integrated  from  b to  a.  The 
variables  and  their  first  derivatives  in  the  system  of 
second  order  equations  for  which  initial  values  are  given 
are  expressed  as  linear  functions  of  the  variables  for 


which  initial  values  are  not  given.  These  expressions  form 
a system  of  n first-order  differential  equations.  The 
coefficients  in  these  expressions  are  the  field  functions. 
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Through  a series  of  differentiations  and  substitutions  of  the 
n first-order  equations,  the  n(n+l)  first-order  field 
equations  are  formed.  These  equations  are  integrated  from 
a to  b to  determine  the  field  functions.  Then  the  n first- 
order  equations  are  integrated  from  b to  a to  determine  the 
solution  to  the  original  problem. 

Although  the  computed  solutions  are  forced  to 
satisfy  the  boundary  conditions  exactly,  it  has  been  shown 
that  with  certain  boundary  conditions  the  solutions  for  the 
internal  points  could  be  greatly  in  error  (Zarghamee,  1970)  . 
This  would  be  due  to  a growing  solution  for  the  field 
equations  when  integrated  from  a to  b.  In  addition,  the 
integration  from  b to  a would  yield  erroneous  results, 
because  the  field  functions  at  b would  be  in  error  and  the 
integration  from  b to  a would  also  give  rise  to  a growing 
solution. 

In  the  invariant  imbedding  method  (Scott,  1973), 
which  is  similar  to  the  field  method,  an  nth-order  differen- 
tial equation,  with  boundaries  at  a and  b,  is  first  replaced 
with  a set  of  first-order  differential  equations.  From  this 
set  of  n first-order  equations  and  with  the  use  of  a modified 
Riccati  transformation,  two  systems  of  equations  involving 
the  invariant  imbedding  fucntions  are  derived.  One  system 
contains  first-order  equations  in  terms  of  invariant  imbedding 
functions  only,  and  is  integrated  from  a to  b to  determine 
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1 the  values  of  these  functions.  The  other  system  contains 

only  algebraic  equations  in  terms  of  the  invariant  im- 
j bedding  functions  and  the  dependent  variables  in  the  original 

system  of  n first-order  equations.  This  system  is  used  to 
determine  the  values  of  the  original  dependent  variables 
between  a and  b.  Here,  however,  is  where  a major  problem 
can  develop.  Expressions  for  these  variables  involve 
taking  differences  of  invariant  imbedding  functions.  In 
certain  cases  these  differences  can  be  very  small  numbers. 

On  the  other  hand,  the  invariant  imbedding  functions  can 
be  very  large  numbers;  therefore,  if  not  enough  significant 
figures  are  carried  the  solutions  could  be  very  much  in  error. 

1. 3 Organization  of  the  Study 

In  Chapter  2,  the  proposed  method  of  analysis  is  developed. 
Section  2.1  gives  a general  description  of  the  process  of  re- 
duction of  order  and  its  use  to  determine  solutions  of  two- 
point  boundary  value  problems  when  growing  solutions  are 
present.  Sections  2.2,  2.3  and  2.4  deal  with  the  explanation 
( of  how  numerical  reduction  of  order  of  differential  equations 

can  be  effected.  First,  a single  second-order  differential 
equation  is  used  to  illustrate  the  reduction  of  order  process. 
The  method  is  then  extended  to  a general  differential  equation 
of  order  n and  finally  applied  to  a system  of  n first-order 
differential  equations. 

i I 


; 
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It  is  shown  in  section  2.4  how  to  use  initial 
value  solutions  and  the  reduction  of  order  procedure  to 
eliminate  difficulties  presented  by  the  growing  solutions. 

The  treatment  of  non-homogeneous  cases  is  also  discussed. 

Chapter  3 presents  an  example  of  a two-point 
boundary  value  problem  which  contains  four  growing  solutions, 
two  in  each  direction.  The  reduction  of  order  procedure 
here  is  used  in  the  determination  of  four  independent  growing 
solutions.  Once  their  values  are  found,  they  are  then 
linearly  combined  to  find  the  correct  solution  which  satis- 
fies the  boundary  conditions. 

In  Chapter  4,  the  natural  frequencies  and  corre- 
sponding mode  shapes  of  spherical  shells  are  determined  by 
the  Holzer  method  (Holzer,  1921)  and  by  a modification  of 
the  Newton-Raphson  method  developed  by  Robinson  and  Harris 
(Robinson  and  Harris,  1971).  The  application  of  both 
methods  involve  the  solution  of  two-point  boundary  value 
problems.  These  problems  are  solved  using  an  initial  value 
procedure  and  the  reduction  of  order  method. 

The  results  of  the  applications  of  the  method  of 
analysis  are  presented  and  discussed  in  Chapter  5 along  with 
results  obtained  from  other  methods.  In  particular,  com- 
parisons are  made  between  the  proposed  method  and  Zarghamee's 
method  for  determining  natural  frequencies  of  spherical  shells. 

In  Chapter  6,  general  conclusions  are  drawn  from 
the  various  applications  of  the  proposed  method  and  possible 
extensions  of  the  procedure  are  discussed. 
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1 . 4 Nomenclature 


Each  symbol  is  explained  when  it  is  introduced  in  the 
text.  The  following  is  presented  for  the  convenience  of  the 


reader. 


Shell  radius.  Chapter  1 


Boundary  conditions.  Chapter  1 


Weighting  coefficient.  Chapter  3 


Linear  differential  operator. 
Chapter  4 


Coefficients  in  general  system  of 
n first-order  differential  equa- 
tions, Chapter  3 


Boundary  condition 


Linear  differential  operator. 
Chapter  4 


Constant  in  differential  equations 
of  reduced  order 


Coefficients  in  fourth-order  system 
of  differential  equations.  Chapter  3 


■ 
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Eh' 


12  (l-v“) 
B and  C 


Y~  / Chapter  4 and  Appendices 


Matrix  of  the  values  of  the  displace- 
ments and  stress  resultants.  Appendix 
0 

Modulus  of  elasticity 

Eh 

(1-v2) 

Right  side  of  general  system  of  n 
first-order  differential  equations 

Shear  force  on  two  adjacent  shell 
segments,  see  Fig.  8 


Sin  <J>  Fn  d6 


Shell  thickness, Chapter  4 and 
Appendices  B and  C 

interval  size.  Chapter  2 


Coefficients  in  fourth-order 


kl'k2 


L 


L 


Li 


n 


n 


M0'M0'M4.8 


differential  equations.  Chapters 
3 and  5 

Interval  length.  Chapter  5 

Differential  operator.  Section  2.2 

Linear  operator.  Section  2.3 

Order  of  system  of  differential  equa- 
tions, Chapters  2 and  3 

Modal  wave  number.  Chapter  4 and 
Appendices  B and  C 

Moment  stress  resultants 


M<pn'M0n'M<t>6n 


Moment  stress  resultants  corresponding 
to  the  value  n,-  Eqs.  (4.1) 
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VNe 


N<J>n'N0n,N4>8n 


P 

P1'P2 


Normal  stress  resultants 

Normal  stress  resultants  cor- 
responding to  the  value  n, 

Eqs.  (4.1) 

Differential  operator.  Chapter  2 

Particular  solutions  to  fourth- 
order  differential  equations 


p(x) 


Right  side  of  fourth-order  differen- 
tial equations.  Chapter  5 


q,r,s,t 


qi'ri'Si,t:i 


ql'rl'Sl' fcl' 
q2'r2'S2'^2 


Dependent  variables  in  the  fourth- 
order  system  of  four  first-order- 
differential  equations.  Chapters 
3 and  5 

Homogeneous  solution  to  the  fourth- 
order  system  of  differential  equa- 
tions, Chapter  3 

First  and  second  homogeneous  solu- 
tions respectively,  to  the  fourth- 
order  system  of  differential  equa- 
tions, Chapters  3 and  5. 
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qpl'rpl'Spl,tpl 

qp2,rp2,Sp2,tp2 


qTl,rTl'STl#tTl 
qT2 ' rT2 ' ST2 ' fcT2 


q 


ij 


Wqz 


First  and  second  particular  solu- 
tions repectively,  to  the  fourth- 
order  system  of  differential 
equations.  Chapter  3 

First  and  second  complete  solutions 
to  the  fourth-order  system  of  dif- 
ferential equations,  Chapters  3 and  5 

Dependent  variable  in  the  general 
system  of  n first-order  differen- 
tial equations.  Chapter  2 

Solution  to  general  system  of  n 
first-order  differential  equations. 
Chapter  2 

Components  of  the  external  pressure 
on  the  shell 


VQe 


Shear  stress  resultants 


r 


a + z 


R(i) 


Residual  function  in  the  Robinson- 


Harris  method 
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I 

I 

r 

i: 

i 

u,^ 

I 

u,  v,w 


u,v,w 


dw 

ad$ 

Time  coordinate.  Chapter  4 

Dependent  variable  in  the  differ- 
ential equations.  Chapter  2 


Components  of  the  displacement 
vector  of  the  middle  surface  in 
the  meridional,  circumferential, 
and  normal  directions  respectively. 
Chapter  4 and  Apendices  B and  C 

Components  of  the  displacement 
vector  of  any  arbitrary  point  in 
the  shell, not  necessarily  on  the 
middle  surface 


Displacement  function  of  the  middle 
surface  corresponding  to  the  value 
n,  Eqs.  (4.1) 


3u,  3v,  3w 
3<t>  3<J>  3 <f> 


respectively 
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v,  v. 
l 


Dependent  variables  in  the  dif- 
ferential equations  of  reduced 
order.  Chapters  2 and  3 


V, 


Qa  + 


n 


4>  a sm  <t>  M 


4>  e 


M 


N 


<j>e 


(j>e 


w,W 


Dependent  variables  in  general 
systems  of  nth-order  differential 
equations.  Chapters  1,  2,  3 and  5 


w_ 


Particular  solution  of  the  fourth- 
order  system  of  differential  equa- 
tions, Chapter  5 


- 


x,x 


(i) 


Eigenfunctions  in  the  Robinson-Harris 
method.  Section  4.5.2  and  Appendix  D 

Dependent  variable  in  general  system 
of  differential  equations.  Chapter  2 

Dependent  variable  in  shell  differ- 
ential equations  of  reduced  order. 
Appendix  B 
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I 

I 

I 

I 

1 

I 

t 


Thickness  coordinate  measured 
positive  outward  from  the  mid- 
dle surface  of  the  shell 

In  (sin  4>)  , Appendix  C 

Boundary  conditions.  Chapters  2 
and  3 


r 

[ 

t 

i 


Weighting  coefficients.  Chapter  3 

Constant  coefficient  in  fourth- 
order  differential  equation. 
Chapter  1 

When  preceeding  a variable,  it 
indicates  an  incremental  change 
in  that  variable 

Determinate  of  the  coefficients  of 
the  linear,  homogeneous  equations 
of  the  boundary  conditions  at  the 
base  and  apex  of  the  shell 
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V'Yer'Y*e 


e<t>,ee'er 


Shear  strains 


Strains 


! 


V 


Poisson's  ratio 


P 


P.O,T 


a 


<D 


o 

r 


x<pr ' T$r ' x 4>0 


<t> 

yo 


<t>,e 


Mass  density  of  shell  material. 
Chapter  4 and  Appendices  A,  B 
and  C 

Dependent  variables  in  differ*,  tial 
equations  of  reduced  order.  Chapters 
3 and  5 

Stresses 

Shear  stresses 

Opening  angle  of  a spherical  shell 

Geometric  coordinates , colatitude 
and  azimuth  angle,  locating  points 
on  the  middle  surface  of  the  shell 


ii 

;i 

ii 

:! 


- «fc 


I 

I 

I w natural  frequency  (rad. /sec.) 

w Reference  frequency  = ~\  / — ^ E _ 

Vpa^(l-V) 

! 

r 

i 

i 

i 

i 

r 

r 
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2.  DEVELOPMENT  OF  THE  METHOD  OF  ANALYSIS 


2 . 1 General 

Solutions  of  differential  equations  occurring  in 
various  types  of  two-point  boundary  value  problems  quite  often 
contain  rapidly  growing  solutions  which  will  dominate  in  any 
initial  value  solution.  Consequently,  solutions  obtained  by 
straight  forward  numerical  integreation  of  the  differential 
equation,  even  if  using  independent  initial  conditions,  will 
not  be  independent  enough  for  numerical  accuracy.  As  mention- 
ed in  Chapter  1,  if  the  solutions  are  not  linearly  independent, 
they  cannot  be  accurately  combined  to  form  a solution  which 
satisfies  the  boundary  conditions. 

The  method  of  reduction  of  order  of  differential 
equation  uses  one  growing  solution  to  obtain  other  solutions, 
either  growing  or  non-growing,  which  are  truly  independent. 

Thus,  the  solutions  can  be  linearly  combined  to  form  a solu- 
tion which  satisfies  the  boundary  conditions. 

In  Section  2.2,  the  well-known  process  of  reduction 
of  order  of  a differential  equation  will  first  be  illustrated 
using  a homogeneous  second-order  differential  equation.  Then, 
the  method  will  be  applied  to  an  nth-order  differential  equation. 
It  will  be  seen  that  difficulties  are  encountered  when  using  the 
straightforward  method  of  reduction  of  order  if  the  first  solu- 
tion obtained  passes  through  zero.  A modification  to  eliminate 
this  difficulty  will  also  be  explained  in  Section  2.2. 
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The  process  of  reduction  of  order  is  extended  to 
the  general  case  of  an  nth-order  system  of  differential 
equations  in  Section  2.3.  This  new  method  is  applied  to  an  nth- 
order  system  consisting  of  n simultaneous  first-order  homo- 
geneous differential  equations  where  the  first  solution  may 
pass  through  zero.  The  process  of  determining  a non-growing 
particular  solution  of  the  general  nth-order  system  will  be 
explained  in  Section  2.4. 

Section  2.5  sets  forth  the  general  process  of  ob- 
taining the  solution  of  a boundary  value  problem  where  the 
corresponding  initial  value  problem  has  growing  solutions. 

It  will  be  seen  that  for  a system  of  differential  equations 
of  order  n,  n+2  initial  value  problems  must  be  solved,  n/2 
homogeneous  ones  from  each  end  (if  n is  even)  and  a parti- 
cular solution  from  each  end. 

2 . 2 Method  of  Reducing  the  Order  of  a Differential  Equation 

Consider  a linear  ordinary  differential  equation 
L (u)  =0  of  order  n,  with  r independent  solutions  known.  It 
is  well  known  (Ince,  1956)  that  by  using  r known  independent 
solutions,  the  equation  can  be  reduced  in  order  from  n to  n-r 
with  only  n-r  solutions  remaining.  To  clarify  this  procedure 
and  show  the  nature  of  the  solutions,  a second-order  system 
will  be  studied  first. 
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Let  be  the  known  solution  of  the  second-order 
differential  equation 


u"  - u = 0 


(2.1) 


corresponding  to  the  following  set  of  initial  conditions: 


u^O)  = a,  u^(0)  =0 


Now  let 


v 


1 


u 

1 


and 


(2.2) 


where  U2  is  a second,  unknown,  solution  of  equation  (2.1) 
and  is  also  unknown.  Solving  (2.2)  for  u2  yields 


dx 


* 


(2.3) 


x 


The  integral  is  a definite  integral,  say  qJ  vdx.  However, 

in  some  cases  it  is  convenient  to  have  a different  lower  limit 
Therefore,  the  definite  integral  will  be  written  as  if  it  were 
an  indefinite  integral. 
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If  we  differentiate  (2.3),  we  get 


u2  = u^v  + u,  / vdx 


(2.4a) 


A second  differentiation  of  (2.3)  yields; 


u2  = 2u1v2  + ulv2  + U1  /v2dx 


(2.4b) 


Substitution  of  (2.3)  and  (2.4)  into  Eq.  (2.1)  gives; 


2u1V2  + U1V2  + U1  /v2dx  " U1  /V2dX  = 0 


or 


/' 


v2dx(u1  - u)  + 2u-lv2  + uv2  = 0 


(2.5) 


Since  u^  is  a solution  of  Eq.  (2.1),  Eq.  (2.5)  reduces  to; 


2u1v2  + uxv2  = 0 


(2.6) 


Thus,  Eq.  (2.1)  of  second  order  has  been  reduced  to  Eq.  (2.6) 
of  first  order  with  v2  as  the  unknown  dependent  variable,  and 
u^  and  u^  as  known  coefficients.  Equation  (2.6)  can  now  be 
integrated  to  determine  v2*  The  initial  condition  for  the 
integration  is  determined  by  assuming  another,  and  independent. 
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set  of  initial  conditions  for  Eq.  (2.1).  Let  the  second  set 
of  initial  conditions  for  Eq.  (2.1)  be: 

u2(0)  = 0,  u2 ( 0 ) = 6 

We  see  from  Eq.  (2.4a)  that  at  x = 0 

u2(0)  = u1(0)v2(0)  ( fv 2dx  = 0 at  x = 0) 


or 


v2(0) 


u2(0)  = 3 
u-^O)  a 


After  v2  has  been  determined  the  other  solution,  u2,  of  Eq. 
(2.1)  is  obtained  by  use  of  Eq.  (2.3). 

In  general,  the  nth-order  equation  is  not  solvable 
analytically  in  terms  of  familiar  functions.  If  a numerical 
integration  of  the  differential  equation  is  carried  out,  then 
u^  will  automatically  include  a significant  component  of  the 
most  rapidly  growing  solution.  Thus,  when  the  equation  is 
reduced  in  order,  this  rapidly  growing  solution  will  no  long- 
er be  a solution  of  the  equation  of  lower  order. 

The  general  case  of  an  nth-order  differential  equa- 
tion L (u)  = 0 is  treated  as  follows  (Ince,  1956) : 


— 
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Let 


_ d 
vr  dx 


vn  = 


n+1 


_d_ 

dx 


_d_ 

dx 


1 . _d_  / 1 

Vr-1  dxVvr_; 


1 . _d_  / 1 

v . dx Iv  . 
n-I  ' n~4 


n 


d_ 

dx 


d 

v2dx 


d 

d 

v^3x 


(2.7) 


where 


L (u)  = 0 

Since  the  order  of  the  differential  equation  is  n,  only  n 

independent  solutions  exist.  Consequently,  u must  be  some 

linear  combination  of  the  n independent  solutions  u^ , Uj» 

. . . u . It  follows  then,  that  v has  to  be  zero, 
n n+l 


Thus 
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° = ~ 
dx 


_L  A ( 1 

v dx  \ v 
n n-1 


d . u \ 

v?*  vl/ 


(2.8) 


If  the  solutions  u^/  . .,  ur, (r<n)are  known  then,  Eq.  (2.8) 

may  be  written  as 


P(v) 


(2.9) 


where  P is  the  differential  operator  of  order  n-r  and 


v = 


dx  vrdx 


v2dx 


— (2.10) 

v 


solving  Eq.  (2.10)  for  u gives 


f , 

> / 

•* 

1 

u = vlJ  \ 

1 V2/ 

V3  * ' * • KVr 

/ vdx]  .... 

J ' 

dx  ) 

J \ 

1 J 

J 

J 

Difficulties  arise  in  the  above  procedure  if  at 
some  point  any  v..  passes  through  zero.  If  this  happens, 
it  is  obvious  that  the  remaining  n-i  solutions  cannot  be 
calculated  at  that  point  by  numerical  integration  without 

t 

introducing  signularities  in  the  remaining  vg,  v^+1,  vi+2* 

. . . •*  vn. 

For  example,  with  different  initial  conditions,  at  some 
point  where 


u^  = v^  = 0 and  u2  / 0 


4_ 

r 


is  undefined.  Any  attempt  to  calculate  u2  at  this  point  us- 
ing the  above  expression  would  involve  numerical  difficulties 
A new  method,  therefore,  had  to  be  devised  for  handling  this 
particular  situation.  A constant  c (c>0)  can  be  added  to  u^ 
such  that  v^  = u^  + c is  always  greater  than  zero.  This  will 
of  course,  change  the  procedure  for  determining  v2,  but  will 
not  change  the  final  complete  solution  of  Eq.  (2.1).  Going 
back  to  the  previous  second  order  example  we  have: 

u"  - u = 0 and  v^  = u^  as  before. 

However,  now 


u2  - (ul  *')''!*  ul/v2dx 

and 

m ii  r i / \ i 

u2  = ui  / v2dx  + 2u^v2  + ^ u 1 + cj  v2  (2.12) 

where  v2  is  different  from  v2- 


(ui + ■=)  / 


v2dx 


(2.11) 


Substituting  Eqs.  (2.11)  and  (2.12)  into  Eq.  (2.1)  and 
simplifying,  we  get 

2u^v^  + (u1  + c)  v2  + c Jv2 dx  = 0 (2.13) 


Thus,  again  we  have  reduced  Eq.  (2.1)  to  a first  order  dif- 
ferential equation.  Here  we  see  that  if  c = 0,  Eq.  (2.13) 
is  the  same  as  Eq.  (2.6).  With  the  occurrence  of  the  term 
cfv^dx  in  Eq.  (2.13)  one  would  tend  to  think  this  equation 
is  of  second  order.  However,  because  J v^dx  is  a definite 
integral,  Eq.  (2.13)  is  of  first  order;  that  is,  it  has 
a one  parameter  family  of  solutions. 

The  initial  conditions  for  Eq.  (2.13)  are  dictated 
by  the  initial  conditions  for  Eq.  (2.1)  as  in  the  earlier  ex- 
ample. The  value  of  J vdx  at  the  initial  boundary  is  zero  as 
before.  The  value  of  v2(0)  is  determined  by  substituting  into 
Eq.  (2.13)  the  independent  sets  of  initial  conditions  for  Eq. 
(2.1),  and  in  this  case  is 


v2(°) 


3 


(a+c) 


(2.14) 


This  is  the  only  initial  condition  possible  for  Eq.  (2.13) 
given  the  two  sets  of  initial  conditions  for  Eq.  (2.1). 
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2 . 3 Reduction  of  Order  of  General  Systems  of  Homogeneous. 
Ordinary,  First-Order  Differential  Equations 

i ( 

A system  of  first  order  differential  equations  is 
generally  easier  to  handle  than  one  differential  equation  of 
higher  order  (Goldberg,  1967)  . The  conversion  of  a differen- 
tial equation  of  higher  order  to  a system  of  first  order  dif- 
ferential equations  merely  involves  the  replacement  of  the  de- 
pendent variable  and  its  derivatives  by  new  variables. 

Let 

(n)  (n-1)  " ' 

y'  + A y'n  + . . . + A y + A y + A y = 0 (2.15) 

i n-2  n-1 , n 

be  an  nth-order  differential  equation  with  the  following 
boundary  conditions : 


y (0) 

= °L 

y(n-2)(L)  - an/2+l 

y'(0) 

= °2 
• 

• 

y(n/2+i)  (L)  . %/2+2 

1 

y((n-2./2)(0) 

• 

• 

• 

= ai\/2 

y (n-i)  (L)  = ^ 

(2.16) 

Here,  n is  assummed  to  be  even,  because  most  engineering  pro- 
blems involve  differential  equations  of  even  order. 


However,  n could  be  odd  in  which  case  (n-l)/2  conditions  would 
be  specified  at  one  boundary  and  (n+l)/2  would  be  specified  at 
the  other. 
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Define  new  variables  as  follows: 


ui = y 


u2  = ui  = y 


U3  = U2  " y 


' * (n-1) 

u = u , = y 
n n-1  1 


• (n) 

un  * y 


(2.17) 


If  Eqs . (2.17)  are  substituted  into  Eq.  (2.15),  we  get  one 
first  order  differential  equation. 


u + A.u  n + . . . + A ,u_  + A u,  = 0 
n ln-1  n-12  nl 


The  other  n-1  first  order  equations  are  taken  from  Eqs. (2.17): 


U1  - u2  - 0 


U2  - u3  * 0 


Ui  ' Ui-1  " ° 


u„  , - u„  = o 
n-1  n 


(2.18) 


. - ■ ■ • v-‘  •>- 

V ‘ -u  ' 


• v^»r  •• 
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The  boundary  conditions  for  the  first-order  system  of  equa- 
tions would  be  the  same  as  the  boundary  conditions  for  the 
nth-order  single  differential  equation:  that  is: 

V0>  - y<°>  - “i  u„/2«(l|  ' yln/2)(U  - V2+1 

I 

u2(0)  = y (0)  = a2  


Un/2(01  = y1"72”1’ (0)  = cn/2  un(L)  - Y(n'1) (L)  - cn  (2.19) 

The  general  system  of  linear  ordinary  first  order 
differential  equations  should  be  written  in  the  following 
specified  form,  solved  for  the  derivatives,  in  order  that 
numerical  integration  techniques  be  readily  applicable: 


ql  * Ll(ql'  q2 V ' 0 


q2  + L2(ql'  q2'  * * *'  qn}  “ 0 


qn  + W q2 qn>  = 0 


(2.  20) 


where  L. , L_,  . . .,  L are  linear  in  the  q's. 
1 z n 


32 


The  boundary  conditions  are: 

qx  ( 0)  = 

qn/2+l (L)  6n/2+l 

q2(o>  = e2 
• 

• 

• * 

qn/2+2(L)  = en/2+2 

qn/2<0)  ' 8„/2 

q (l)  = en 

n n 

(2. 

Assume  that  q , q21#  . . . . 

, q^  is  a growing  solution 

mined  by  straight  forward  numerical  integration  of  Eq.  (2.19) 
starting  with  the  following  set  of  assumed  initial  conditions: 

qu<01  = Cj 

q21(°l  ; o2 

qnl(°)  = % (2.22) 


.'n 


T 


Let 


where 

mined; 
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q12  = (qll 

q22  = q21 J 

q32  = q31 f 

qn2  qnl f 

q12 ' q22 ' * 
that  is. 


h 


+ c)  /v.dx 


v^dx  + v 2 


v^dx  + Vj 


v,dx  + v 
1 n 


*n2 


( q 1 1 (x)  + c)  > 0 
and  (c  > 0) 

o < x £ = L 

(2.23) 

is  a second  solution  to  be  deter- 


q12  + Ll(q12,  q22 qn2>  = 0 

f 

q22  + L2  (qi2,  q22 , * * ‘ *'  qn2)  = 0 


qn2  + Ln  (q12'  q22"  * * *'  qn2)  " 0 


Substitution  of  Eqs.  (2.23)  and  their  derivatives  into  Eqs 
(2.24)  yields 
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(qll  * c)  v!  * qn/vldx  + L1  [(qll  + °*  /vldx' 

q21  fv ldx  + v2 q„l/vldx  * v„]  = 0 

q21Vl  * q21 y ldx  * v2  * L2  (qll  + c>  ldx- 

q21  JVldx  + V2 q„l/vldx  + VJ  * 0 


qnlVl  + V /vldx  + \ + Ln  [<qn  + c>  /* ldx' 

q21  Jv ldx  + v2 qnl/vldx  * vnj  - 0 

Rearranging  terms,  we  get 


,<511'  q21'  ‘ ' 


(qll  + C,V1  + ?11+  Vqll'  q21' 

+ 1^  c v1dx,  v2,  v3, 

LV1  + v2  + [ q21  + Vqll'  q21' 
+ L*  [°  VldX'  ^ 


L2 (qll'  q21'  * ' 

f V2  / V3  , • • 


• •'  qnl>]  /V 

• •'  v„  J * 0 

• " qnl>]  /v 


' " qnl) 
• •'  vn] 


qnlV 


• r • 

1 + Vn  + [ qnl  + Ln 
+ L |c  v.d 

n l 1 


(qll'  q21 ' 


*'  qnl} 


] /-, 


■jdx,  v2,  v3. 


However, 
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<qn  + C)vi  + 4 

[°J 

fvid*' 

V2 ' 

V3 ' ' 

. . . , V 

n 

= 0 

V2  + Wl  + L2 

b 

fv.dx, 

v2 ' 

V3 ' • 

. . . , V 

n 

= 0 

1 

v + q ,v,  + L 

'.j 

• • • 

r 

v,  dx. 

v_ , 

v_ , . 

. . . , V 

= 0 

n ^nl  1 n 

J 

1 

2' 

3 

n 

Therefore, 

<qu  + c)Vl  + 4 

LCJ 

r 

Vjdx, 

V2 ' 

V3 ' * 

. . . , V 

n 

= 0 

V2  + q21Vl  + L2 

,cj 

fv^x, 

V 2 ' 

V3'  * 

. . . , V 

n 

= 0 

, r 

r 

v + q -V  + L j 

n ^nl  1 n |_ 

cJVfdx,  v2,  v3 vn 

I 


Hence,  the  system  of  dif f erential  Eqs . (2.20)  of  order  n has 

been  reduced  to  the  system  of  differential  Eqs.  (2.25)  or  order 

n-1.  The  initial  conditions  for  the  equations  of  reduced 

order,  as  previously  mentioned  in  Chapter  1,  are  dictated  by 

the  chosen  independent  sets  of  initial  conditions  for  Eqs. 

(2.20).  Equations  (2.25)  can  now  be  numerically  integrated 

to  find  the  variables  v. , . . . . , v . Once  their  values 

i n 

are  obtained  another  solution  of  the  system  (2.20)  can  be 

* 

determined  by  use  of  Eqs.  (2.23)  . 


Further  explanation  of  the  numerical  procedure  used  to  calcul- 
ate the  second  solution  will  be  given  in  Chapter  3. 


1* 
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If  system  (2.20)  contains  only  one  rapidly  growing 
solution,  it  will  be  ^l'  *321'  • • ♦ •»  ^ni*  Other  solutions 
of  (2.20)  will  not  grow  rapidly.  On  the  other  hand,  if  the 
system  (2.20)  contains  more  than  one  rapidly  growing  solution, 
then  the  second  solution  q^/  ^22'  • • • ••  *3n2  he  a 

growing  solution  also,  but  independent  of  the  first  and  more 
slowly  growing  than  the  first.  This  second  solution  can  be 
used  to  reduce  the  order  of  Eq.  (2.25)  and  a third  solution  of 
Eq.  (2.20)  can  be  determined.  The  process  of  reducing  the 
order  should  continue  until  all  growing  solutions  have  been 
obtained. 

2 . 4 Procedure  for  Solving  nth-Order  Systems  Using  Reduction 
of  Order 

In  the  second-order  case,  independent  initial 
value  solutions  for  the  homogeneous  equation  can  be  deter- 
mined without  reduction  of  order.  Numerical  integration 
of  the  homogeneous  second-order  differential  equation  will 
yield  independent  solutions,  if  one  solution  is  obtained  by 
integrating  from  0 to  L and  the  other  from  L to  0.  An  in- 
dependent slowly  varying  particular  solution  cannot,  however, 
be  obtained  without  reducing  the  order  of  the  equation.  In 
general,  for  a system  of  differential  equations  of  order  n 
containing  rapidly  growing  solutions,  a reduction  of  order  is 
necessary  to  obtain  independent  homogeneous  solutions  as  well 
as  an  independent  particular  solution. 
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A system  of  differential  equations  of  order  n, 
where  n is  even,  could  possibly  have  n growing  solutions. 
Commonly,  n/2  solutions  grow  from  0 to  L and  n/2 
solutions  grow  from  L to  0 . Unless  it  is  known  in  ad- 
vance how  many  growing  solutions  are  present,  it  is  neces- 
sary to  keep  reducing  the  order  of  the  system  of  differen- 
tial equations  to  find  each  new  solution. 

The  procedure  for  determining  the  initial  value 
solutions  for  a general  system  of  linear  ordinary  differen- 
tial equations  of  even  order  n containing  n growing  solutions 
is  as  follows: 

1.  Integrate  the  original  nth-order  system 
of  differential  equations  from  0 to  L to 
obtain  the  fastest  growing  solution  from 

0 to  L,  u^. 

2.  Use  u^  to  reduce  the  order  of  the  equa- 
tion to  n-1  and  integrate  the  reduced  order 
equations  from  0 to  L. 

3.  Use  the  solution  of  the  n-1  order  system 
of  equations, v^,  and  the  first  solution  of 
the  nth  order  differential  equation,  u1 , to 
determine  another  solution,  u2,  of  the 
original  nth  order  system  as  previously  des- 
cribed . 
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4.  Test  the  second  solution  to  the  nth  order 
system  to  see  if  it  grows  rapidly.  If  it 
does  not  grow  rapidly,  then  no  further  re- 
ductions of  order  are  necessary  to  determine 
the  remaining  independent  solutions.  If  it 
does  grow  rapidly,  then  another  reduction  of 
order  is  necessary  to  obtain  another  indepen- 
dent solution  of  the  nth-order  system. 

5.  Assuming  further  reductions  of  order  are 
necessary,  use  solution  v^,  already  obtained, 
to  reduce  the  n-1  order  system  to  an  n-2  order 
system.  Solve  this  system  and  use  its  solu- 
tion, V2»  along  with  solutions  u^  and  v^  to 
obtain  solution  number  3,  u^,  of  the  nth  order 
system. 

6.  In  general,  the  n-th  solution,  un_i*  of 
the  original  nth-order  system  is  determined 
using  solution  u^  of  the  nth-order  system  and 
solutions  v^,  . . . .,  vn+i_i  of  the  systems 
of  reduced  order. 

7.  If,  at  some  point,  a solution  is  obtained 
that  does  not  grow  rapidly,  further  reductions 
are  not  necessary  to  determine  the  remaining 
initial  value  solutions  from  0 to  L.  The  re- 
maining solutions  are  then  obtained  by 
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integrating  the  latest  reduced  order  equation 
from  0 to  L a number  of  times  starting  with 
independent  sets  of  initial  conditions. 


The  initial  value  solutions  which  grow  from 
L to  0 are  obtained  in  the  same  manner  as  the 
solutions  which  grow  from  0 to  L,  however,  the 
direction  of  integration  is  reversed. 


2.5  Non-growing  Particular  Solutions 


Up  until  now  we  have  only  discussed  the  numerical 
determination  of  independent  homogeneous  solutions.  Many  two- 
point  boundary  value  problems  involve  differential  equations 
which  are  non-homogeneous . It  then  becomes  necessary  to 
determine  an  independent  particular  solution  as  well  as  a 
full  set  of  independent  homogeneous  solutions. 

A particular  solution  can  be  obtained  using  the 
same  procedure  as  for  determining  solutions  (Ji2'  • • • •# 
q^n^2*  However,  in  this  case  the  reduced  order  equations 
and  the  original  first  order  system  of  differential  equations 
will,  of  course,  be  non-homogeneous.  Equations  (2.20)  become 


q.  + L.(q. , q-,  . . . .,  q ) ■ f,(x) 
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If  the  right  side  of  Eq.  (2.26)  is  carried  along  in  the  re- 
duction process,  we  find 


(q,!  + c) 


V1  * hllcfvl ' v2 V - £l' 


2 * q21Vl  + L2  <c/v  V2' 


. V ) = f_(x) 

n 2 


vn  + qnlvl  + VCJV  V 


" Vn)  = fn(x)  (2.27 ) 


where  q^,  . . . q ^ is  the  original  homogeneous  initial 
value  solution  to  Eq.  (2.20).  From  this  point  on,  the  solu- 
tion is  computed  exactly  the  same  way  as  the  other  solutions 
but  using,  for  convenience,  homogeneous  initial  conditions. 

The  reductions  are  done  i times  until  the  reduced 
homogeneous  equation  of  order  n-i  does  not  contain  a grow- 
ing solution.  This  n-ith  order  equation  will  be  used  along 
with  the  homogeneous  solutions  to  determine  a non-growing 
particular  solution,  Up^.  To  assure  accuracy  of  the  total 
solution  to  the  nth-order  equation  throughout  the  interval 
0 to  L,  it  is  convenient  to  find  a second  particular  solu- 
tion starting  from  L and  integrating  to  zero.  This  particular 
solution  will  involve  the  homogeneous  solutions  which  were 
found  by  integrating  from  L to  0.  The  two  particular  solu- 
tions are  used,  along  with  the  homogeneous  solutions,  to  deter- 
mine two  solutions  to  the  original  nth-order  differential 
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equation.  The  particular  solution  which  was  found  by  inte- 
grating from  0 to  L,  when  combined  with  the  homogeneous  sol- 
utions, will  yield  a solution  which  is  most  accurate  in  the 
first  half  of  the  interval.  The  second  solution  which  con- 
tains a particular  solution  found  by  integrating  from  L to 
0,  will.be  most  accurate  in  the  second  half  of  the  interval. 
The  two  solutions  will  match  near  the  middle  of  the  interval. 
Thus,  the  total  solution  to  the  nth-order  system  will  be  a 
combination  of  the  first  half  of  one  solution  and  the  second 
half  of  the  other  solution. 

2 . 6 Numerical  Integration  Procedure 

The  specific  numerical  integration  method  used  was 
the  trapezoidal  rule  in  a predictor-corrector  type  procedure. 

A procedure  of  this  form  was  chosen  over  other  methods  such 
as  Runge-Kutta,  because  it  is  easy  to  handle  for  certain  cases 
and  permits  the  changing  of  the  interval  size  during  the  inte- 
gration. Variation  of  the  interval  size  is  done  to  correspond 
to  changes  in  the  gradient  of  the  function.  The  definite  in- 
tegrals in  equations  of  the  form  of  system  (2.23)  were  deter- 
mined by  quadratures  using  Simpson's  Rule. 

2.6.1  Method  of  Integration 

There  are  many  techniques  which  can  be  used  for 
numerical  integration  such  as  Runge-Kutta  and  various  predic- 
tor-corrector methods.  For  the  applications  given  herein,  a 
predictor-corrector  method  will  be  used.  The  algorithm  for 
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the  method  of  integration  is  as  follows: 

1.  The  assumed  and  given  initial  conditions 
are  substituted  into  the  differential  equa- 
tions at  point  1 (x=0)  to  calculate  the 
value  of  the  first  derivative. 

2.  The  value  of  the  first  derivative  at  point 
1 becomes  the  first  approximation  of  the  first 
derivative  at  point  2.  The  value  of  the  vari- 
able will  be  determined  by  numerical  integra- 
tion using  the  trapezoidal  rule;  that  is, 

= qi_1  + 1/2  h^  + qi_1)  . 

3.  Values  of  the  variable  at  point  2 will  be 
substituted  back  into  the  differential  equations 
to  calculate  new  values  for  the  first  derivatives. 

These  new  values  will  be  compared  to  the  old  values. 

If  the  two  values  are  within  a certain  tolerance, 

we  proceed  to  point  3.  If  they  are  not,  then  the 

. 

new  value  of  the  first  derivatives  will  be  used  in 
the  trapezoidal  expression  to  compute  another  value 
of  the  variable  at  point  2. 

4.  The  procedure  is  repeated  until  a certain  toler- 
ance has  been  obtained  between  the  old  and  new  esti- 
mates of  the  first  derivatives. 

0 
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5.  If  it  takes  more  than  six  iterations 
for  the  convergence,  the  size  of  the  in- 
terval is  decreased  and  the  procedure  is 
repeated  at  the  same  point.  If  the  con- 
vergence takes  less  than  three  interations, 
the  size  of  the  interval  is  increased  for 
the  next  interval 
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where  t < t < t • Because  the  error  of  one  interval  is  of 
° — — s 

the  order  h , the  error  over  the  total  range  is  of  the  order 
2 

h since  the  number  of  intervals  equals  T/h  (Thomas  1961)  . 
The  error  in  Simpon's  rule  can  be  expressed  as 


h/3  (x 


s-1 


+ 4_  + x 


s+1 


r s+1 

(x(t 

s-1 


) ) dx 


5 4 

h d x (t) 
90  dt4 


Over  an  extended  range  T the  error  would  be 


h/3(x  + 4x,  + 2x_  + 4x,  + 2x.  + 

o 1 2 3 4 


+ 2 


X2s-2 


f T 25 

+ 4x2 s-1  + X2s>  V X(t 


) dx 


4 4 

Th  a x ( T ) 
180  dt4 


(At  = h = T/2s  and  tQ  < t < t2S) 

Here  the  error  over  one  interval  is  of  the  order  h5  and  the 

4 

error  over  the  entire  range  is  of  the  order  h . 

In  addition  to  the  discretization  error  described 
above,  round-off  error  is  also  present.  Round-off  error  will 
increase  as  the  number  of  intervals  increase.  The  interval 
size  selector  is  used  in  the  procedure  to  save  execution  time 
and  space  by  keeping  the  interval  size  at  a maximum  length 
and  the  number  of  intervals  to  a minimum. 
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3.  AN  EXAMPLE  OF  MULTIPLE  EDGE  EFFECTS 


3. 1 General 

In  this  chapter  the  method  of  reducing  the  order 
of  a differential  equation  will  be  applied  to  a system  of 
differential  equations  containing  four  growing  solutions, 
two  in  each  direction.  A fourth-order  example  is  given 
to  test  the  ability  of  the  procedure  to  separate  growing 
solutions.  The  fourth-order  equation  will  first  be  re- 
written as  a system  of  four  first-order  equations,  as  in 
Chapter  2,  before  the  procedure  begins. 


3. 2 Fourth-Order  System  With  Two  Growing  Solutions  in 
Each  Direction 

It  is  easy  to  write  down  a single  fourth-order 


system  with  two  growing  solutions  from  each  edge.  For  in- 


stance, a general  solution  of  the  homogeneous  equation  having 
the  desired  fcrm  could  be  taken  as 


C^e^ix  + C2e^2X  + C^e'^l*  + C^e~^2X  (3.1) 


It  is  obvious  that  e^lx  and  e^2x  will  grow  from  0 to  L 
and  e ^lx  and  e ^2X  will  grow  from  L to  0.  The  fourth- 
order  differential  equation  which  has  Eq.  (3.1)  as  the  homo- 


geneous part  of  its  solution  is 
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where 


k^  > k2  and  k2L  >>  1 

If  ^ > k2  then  k^  represents  the  more  rapdily  growing 
edge  effect.  Also,  k2L  » 1 in  order  that  the  edge  effects 
be  pronounced.  Expanding  the  above  equation  yields 

4w-  (k^  + k^)~-  w + Alvi  = f (x) 
dx4  1 2 dx* 

or 


wlv  - (k^  + k2)w"  + k^w  = f (x)  (3.2) 

The  fourth-order  equation  (3.2)  can  be  re-written 
as  a system  of  four  first-order  equations  by  substituting  new 
variables  for  the  dependent  variable  and  its  derivatives. 

Let 


q = w 


r = q ' = w ' 


s = r ' = w ' ' 


t = s'  = w* ’ • 


f = w*v 


(3.3) 
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Substitute  Eq.  (3.3)  into  Eq.  (3.2)  to  find  one  of  the  four 
first-order  equations 

. 

t'  " (kl  + k2)s  + klk2q  = f(x)  (3.4a) 

The  other  three  first-order  equations  are 

q'  - r = 0 (3.4b) 

r'  - s = 0 (3.4c) 


s’  - t = 0 (3. 4d) 

3. 3 Reduction  of  Order 

When  the  system  of  four  first-order  differential 
equations,  Eqs.  (3.4),  is  numerically  integrated  from  0 to  L 
the  fastest  growing  solution  in  this  direction  will  dominate 
when  x is  at  all  large.  Let  this  solution  be  q^,  r^,  s , 
Another  solution  can  be  defined  by 

c > 0* 


f 
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i | 


= S1  / 


vdx  + a 


vdx  + T 


(3.5) 


(p,  a and  t are  functions  of  x) 


where 


(q1  + c)  > 0 

Substituting  Eq.  (3.5)  into  Eq.  (3.4)  and  simplifying  we 
find  the  following  thrid-order  system  of  differential  equa- 
tions : 

p'  - (k^  + k^)u  + k^k^c  Jvdx  + tv  = f(x) 
a'  - t + sv  = 0 

P ' - a + rv  = 0 

v (q  + c)  - P = 0 (3.6) 

These  equations  in  turn  can  be  solved  for  P,  cr  , x,  and  v. 
The  second  solution  q2,  r2,  t 2 described  by  Eq.  (3.5) 

will  be  forced  to  be  independent  of  the  first  solution  and 
more  slowly  growing. 


1 

• . 


| 


I 

I 


T 
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In  the  numerical  integration  of  the  original 
nth-order  differential  equation,  it  is  found  that  the 
solution  q^ , r , s^,  t contains  significant  portions 
of  all  the  growing  and  non-growing  homogeneous  solutions 
when  x is  small.  However,  when  x is  large  the  most  rapidly 
growing  solution  that  might  be  contained  in  q2,  r2,  s2,  t 2» 
is  forced  to  equal  zero  at  1.5L,  that  is: 

1.5L 

q2(1.5L)  = 0 = (q^l.SL)  + c)  j vdx  + 0^(1. 5L) 

O 

Thus,  a second,  independent,  solution  is  computed  as  follows 


q2  • (q 


i* c)  r 


vdx  + 


Vl 


r2  = rl*  J vdx  + + P 

s 2 = s^(  J vdx  + ot^)  + o 
t2  = t^(  Jvdx  + c^)  + x 


(3.7) 


By  using  the  above  procedure,  but  integrating  from 
L to  0,  the  fastest  growing  solution  from  L to  0,  q^,  r^, 

s^#  t^,  and  the  slower  growing  solution  from  L to  0,  q^, 

r4,  s4,  t4,  can  be  computed. 

The  reason  that  this  procedure  works  can  be  under- 
stood readily  if  it  is  remembered  that  the  desired  second 
solution  grows  much  more  slowly  than  the  first.  By 
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eliminating  the  first  solution  well  beyond  x = L,  the 
proportion  of  the  first  solution  present  in  the  second 
is  reduced  almost  to  zero  over  the  entire  range  0 to  L. 
Experience  indicates  that  the  combination  at  1.5L  works 
effectively  without  increasing  computational  effort  greatly. 

After  all  four  independent  homogeneous  solutions 
are  found,  the  particular  solution  is  calculated  using  the 
non-homogeneous  third-order  system  (3.6).  To  guarantee 
accuracy  of  the  total  solution  at  both  edges,  two  parti- 
cular solutions  are  determined  and  used  to  calculate  two 
solutions  to  the  nth-order  system. 

One  particular  solution  is  obtained  by  first  in- 
tegrating Eqs.  (3.6)  from  0 to  L and  replacing  q,  r,  s, 
t with  q1#  r^,  s^,  t.^  to  determine  p2,  °2'  T2'  v2‘ 

Then,  by  using  the  following  system,  q , r^  , s , t 

P1  P1 


(3.8) 


pl  pl 


is  determined; 


qpl " 

(ql 

+ c ) Jv 

2dx 

t «2qx  + 

a3q2 

r*r 

rl(. 

Jv2dx  + 

a2) 

+ P2  + 

V2 

s = 
pl 

sl( 

/v2 dx  * 

a2) 

+ a 2 + 

a3S2 

t 

pl 

V 

J v2dx  + 

a2) 

+ T2  + 

a3t2 

I 


where  a 2 and  are  determined  from  the  following  two 


equations : 


(q^ ( L)  + c) J v2dx  + a2  q1(L)  + a3q2(L)  = 0 


:l(W(/ 


v_dx  + a_)  + p _ ( L)  + a,r  (L)  * 0 (3.9) 


The  use  of  Eqs.  (3.9)  guarantees  that  all  rapidly  growing 

solutions  are  eliminated  from  q , r , s , t 

P1  P1  P1  pl  ’ 

A second  particular  solution  is  obtained  using  the 
above  procedure,  but  reversing  the  direction  of  integration 
and  replacing  q,  r,  s,  t in  Eqs.  (3.6)  with  q3,  r3,  s3. 


t3*  The  second  particular  solution  then  is 


(q3  + c)  / v3dx  + V3!  + a5q4 


' r3(/' 


v3dx  + a4)  + P3  + a5r4 


s3(  jv3dx  + a4)  + c3  + a5s4 


V Jv3dx  + «4)  + T3  + °5t4 


where  a4  and  a,.  are  determined  as  follows; 


(3.10) 
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r° 

< q 3 ( 0 ) + c)  / v3  dx  + a4q3(0)  + c^q  4 ( 0 ) = 0 


r3  ( 0)  ( J'vjdx  + a4)  + P3(°)  + a5r4  ( °)  = 0 


The  particular  solution  q , r , s_  , t will  be  most 

P2  ?2  P2  p2 

accurate  when  x is  large. 

After  all  six  solutions  have  been  found,  the  four 
homogeneous  solutions  are  linearly  combined  and  added  to  each 
particular  solution  to  satisfy  the  four  boundary  conditions. 
Assume  that  the  boundary  conditions  are 


q ( 0)  = b1  s ( 0)  = b3 


r (L)  = b2  t(L)  = b4 


and  the  six  solutions  are 

qi'  ri'  si'  ti 

where  i = 1,  2,  3,  4 are  homogeneous  solutions  and  i = 5,  6 
are  the  particular  solutions. 

At  the  boundaries  the  linear  combination  is  the 
solution  and  must,  therefore,  satisfy  the  boundary  conditions. 
Thus, 
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For 

Solution  1 

For 

Solution 

i 2 

4 

4 

iil 

q±  (0) ai 

+ q5(0) 

" "l 

iil 

(0)ai 

+ q6(0)  = b 

4 

4 

iil 

ri  (L) ai 

+ r5(L) 

= b2 

iil 

ri<L)ai 

+ rg(L)  = b2 

4 

4 

iil 

si(0)ai 

+ s5(0) 

" b3 

iil 

si(0)ai 

+ Sg  ( 0)  = b3 

4 

4 

iil 

t±  (L)  ai 

+ t5(L) 

" b4 

iil 

t.  ( L) a . 
i i 

+ tg(L)  = b4 

where  a and  a are  the  weighting  coefficients.  Equa 

tions  (3.11)  are  solved  for  the  unknown  a.  and  a. . The 

1 l 

solution  is  then  found  by  linearly  combining  all  the  solu' 
tions  with  their  respective  weighting  co-efficients  as 
follows : 


Solution  1 

4 

" iil  Vi  + q5 


4 

Z rja:  + r 

i=l 


4 

i£l  8iai  + 3 5 


4 

ill  fciai  + 


Solution  2 


qT,  " ili  qiai  + q6 


rT2  = iil  riai  + r6 


'T  * iil  siai  + s6 


fcT2  = iil  tiai  + fc6 


(3.11) 


total 
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where  qT  , rT  , sT  , tT  and  qT  * rT  , sT  , tT  are  two 

total  solutions.  The  first  solution  will  be  most  accurate 

when  x is  small  and  the  other  solution  will  be  most  accurate 

when  x is  large.  The  two  solutions  will  be  identical  near 

the  middle  of  the  interval.  Thus,  the  total  solution  to 

the  nth-order  system  will  consist  of  q , r_  , s , t 

T1  T1  T1  T1 

for  small  x and  q , r , s , t for  large  x. 

1 2 x2  l2  l2 

Numerical  results  for  fourth-order  example  and 
other  applications  will  be  given  in  Chapter  5. 
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4.  FREE  VIBRATION  OF  THIN  SPHERICAL  SHELLS 


: 

I 


4.1  Introduction 

The  shells  in  this  study  will  be  limited  to  linear- 
ly elastic  spherical  shell  segments.  Love's  first  approxima- 
tion (Love,  1881)  will  be  used.  This  approximation  is  based  on 
the  following  assumptions: 

(1)  The  shell  is  thin  (h/a  1/20)  . 

(2)  The  deflections  are  small. 

(3)  A condition  of  approximately  plane  stress  exist. 

(4)  The  normal  to  the  middle  surface  remains  normal, 
straight  and  short  enough  that  its  extension  can 
be  taken  as  zero  during  deformation. 

A system  of  eight  first-order  differential  equations 
of  motion  for  a spherical  shell  segement  will  be  put  into  a 
usable  first-order  form.  The  equations- will  be  determined 
from  the  following  fundamental  equations  of  a shell  found  in 
the  appendix: 

(1)  The  equations  of  motion  obtained  from 
the  balancing  of  forces  acting  on  some 
fundamental  element. 

(2)  Strain-displacement  relations  obtained 
from  purely  geometric  considerations. 

(3)  The  relationship  between  stress  and 
strain  (Hooke's  Law). 
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4 . 2 Fundamental  Equations  of  a Shell 

The  equilibrium  equations  (A. 1)  and  (A. 2)  in  the 
appendix  can  be  converted  to  a system  of  eight  first-order 
differential  equations.  For  ease  of  calculation,  we  have 
assumed  that  the  motion  of  the  shell  is  symmetric  motion, 
that  is,  it  does  not  include  torsion.  The  displacements 
and  stress  resultants  in  the  hemispherical  shell  are  either 
symmetric  or  anti-symmetric  about  a meridian  from  which  the 
azimuth  angle  0 is  measured. 

Fourier  series  will  be  used  to  change  the  partial 
differential  equations  to  ordinary  differential  equations. 

The  general  expression  of  a function,  say  f(x)/  in  a Fourier 
series  is 

a 

f(x)  = + a^cosnx  + a2cos2nx  + . . . . 

+ b^sinnx  + b2sin2nx  + . . . . 

The  cosine  terms  of  the  Fourier  series  are  symmetrical  and  the 
sine  terms  of  the  Fourier  series  are  anti-symmetrical  about 
zero.  Hence,  the  displacements  and  stress  resultants  can  be 
expressed  in  terms  of  either  sine  or  cosine  functions  of  the 
Fourier  series;  that  is. 
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N. 


N, 


N, 


9 9 


nIo  un  ^ ft)  cos  n 6 

?,  v (4,t)  sin  n 0 
n=l  n 

nI0  wn  M>»t)  cos  n 0 

N (4>,t)  cos  n 0 
n-0  9 n 

? N M>,t)  cos  n 0 
n-1  g n 

% = n-1  N*0n  (<pft)  sin  n 0 


% = nlo  ®<J>  n (*'t}  cos  n 0 

% = „Ii  (*'t}  Sin  n 6 

M*  = nlo  l\n  (*ft)  COS  n 8 

% = nIo  ^n  (*  't}  COS  n 9 

M*0  = = nIlM0?nsin  n 0 

(4.1) 


The  displacements  u and  w are  even  functions;  therefore,  the 
Fourier  series  for  them  only  involves  the  even  functions 
cosn0.  Conversely,  the  displacement  v is  an  odd  function 
and,  therefore,  its  Fourier  series  only  involves  sinn0. 
Likewise,  the  stress  resultants  are  either  even  or  odd 
functions;  thus,  their  Fourier  series  only  involve  either 
sinn0  or  cosn0. 

Substitution  of  the  Fourier  series  expressions  for 
the  displacements  into  Equations  A. 6 yields  the  following: 


Eh  nvn 

N,  = K—  t W + U . + V (w  + — I T-  + u cot  )] 

♦ a (i_v  ) n n'  n sin  <J>  n * J 


Ph  nvn 

He  “ 77r;77  Iw" + nr?  * COt  * un  + (wn  + V ♦>  1 

a i l” v ) 
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(4.2) 


Using  the  above  equations  and  the  equilibrium  equations,  we 
can  derive  a system  of  eight  first-order  differential  equa- 
tions of  motion  of  the  shell. 


Let 
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dw 
a d<J> 


dw 

d<J> 
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(4.3a) 


* 

For  ease 


in  writing, 


the  subscript 


n will  be  dropped. 


ltd 
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This  becomes  the  first  of  the  eight  first-order  differential 
equations.  Through  a series  of  substitutions  involving  the 
stress  resultant  equations,  the  equilibrium  equations,  de- 
rived in  Appendix  A,  and  the  Fourier  series  expressions,  the 
other  seven  first-order  equations  are  determined.  Substitu- 
tion of  Equation  (4.3a)  into  (4. 2d)  yields: 
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a - u, 
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Solving  the  above  equation  for 

ds  a\  . ( 

% = a d$  - ~d  + v[s  a cot  ♦ 
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Substitution  of  this  expression  into  4.2a  gives 
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ds 

Solving  for  a ^ yields 


ds 

a d<J> 


.(  . v"2.  - - - ( - c°i.t ) 3 + 

v a2  sin2  * > \ * ) aE  D 


which  is  the  second  equation  in  the  eighth-order  system. 
Similarly,  the  remaining  six  first-order  differential 
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equations  can  be  found.  They  are  as  follows: 
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where 


E'  = 


Eh 

1-v" 


D = 


Eh' 


12  (1-v  ) 


k = 1 + 


a"^7 


Tne  terms  and  V0  are  similar  to  the  Kirchhoff  shear  terms 

and  are  derived  in  Appendix  E.  The  terms  X,  Y,  Z,  U and  V 
are  inertia  terms  from  the  force  and  moment  equilibrium 


equations . 


The  first  term 


X , is  expressed  as  follows: 
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Neglecting  higher  order  terms  we  get 
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2 a2 
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Y sm  <J>  = p a sin  <*>  -2— =■ 

3t 


2 a2 

Z sin  <p  = p a'n  sin  <J>  -2-S 

at 


U sin  4>  = 0 (h  ) (which  will  be  taken  as  zero) 


V sin  <p  = 0(h  ) 


(also  taken  as  zero) 


In  order  to  determine  the  natural  frequencies  and 
mode  shapes  of  a shell  the  actual  external  loads  have  to  be 
set  to  zero,  and  a harmonic  motion  assumed.  That  is. 


u ( 4> , t ) 
v(d>,t) 
w(<J>,  t) 


v (<j>) 
w (<j>) 


and  thus. 


N(<t>,t) 
s (<J>,  t) 
Q^(<t>,t) 


V#) 

M,  „ (4>) 


(4.4) 


Substitution  of  Equations  (4.4)  into  Equations  (4.3)  yields 
the  following  system  of  eight  first-order  homogeneous  equa- 
tions : 
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(4.5) 


These  equations  will  be  used  to  find  the  natural  frequencies 
and  mode  shapes  of  the  shell. 

4 . 3 Boundary  Conditions 

The  system  of  eight  first-order  shell  equations  con- 
tain eight  unknowns.  Therefore,  it  is  necessary  to  have  eight 
boundary  conditions  to  determine  a unique  solution,  four  at 
the  base  and  four  at  the  apex.  The  boundary  conditions  at 
the  base  are  generally  quite  easy  to  determine.  They  depend 
on  the  type  of  support  of  the  shell.  Here,  all  the  shells  will 
be  taken  to  have  a fixed  base.  This  means  that  at  the  base 
there  can  be  no  displacements  in  the  meridional,  circumferential, 
or  radial  directions.  Also,  the  rate  of  change  of  the  meridional 
displacement  at  the  base  with  respect  to  the  meridional  angle  <p 
must  be  zero.  These  conditions  hold  true  from  t-initial  to 
t-final.  That  is. 
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u ( <J> » t ) = 0 
v(<J>Q,t)  = 0 
w(<J>o,t)  = 0 

- 0 


The  boundary  conditions  at  the  apex  of  a hemispheri- 
cal shell  are  considerably  more  complex  than  those  at  the  base. 
Although  the  actual  physical  quantities  are  not  singular  at 
the  apex,  a singular  point  exists  in  the  coordinate  system 
chosen.  At  the  apex,  a one-to-one  correspondence  between  a 
point  as  described  by  the  coordinate  system  and  an  actual  point 
on  the  shell  breaks  down.  For  example,  the  points  0^,  r) 

and  0^+A0,  r)  are  the  same  point  if  <^=0  regardless  of 

the  value  of  A0. 

We  know  that  the  apex  does  not  behave  any  differently 
than  any  other  part  of  the  shell.  Therefore,  the  stress  result- 
ants, the  displacements  and  their  derivatives  have  to  be  finite 
at  the  apex.  Examination  of  the  eight  first-order  equations 
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Thus,  we  have  three  conditions  for  n=0,  two  for  n=l  and 
four  conditions  for  n>l.  For  the  first  case,  n=0,  only 
three  conditions  at  the  apex  are  necessary,  because  the 
eighth-order  system  of  equations  reduces  to  sixth-order 
when  n=0.  When  n>l  we  have  tour  conditions  at  the  apex 
needed  to  determine  a solution.  However,  for  n=l  we  are 
lacking  two  necessary  boundary  conditions.  Greenbaum  (Green- 
baum,  1963)  added  two  conditions  to  the  case  where  n=l  by 
considering  the  force-strain  relationship  and  moment-curvature 
relationship  of  the  shell  segment.  The  results  were  as  follows 
for  n=l: 


u'  = 0 


These  two  conditions  with  the  previously  determined  conditions 
constitute  the  four  necessary  boundary  conditions  at  the  apex 
for  n=l. 
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4 . 4 Method  of  Numerical  Integration 

The  system  of  shell  equations  derived  in  Section  4.2 
is  of  eighth-order  when  n>l  and  sixth-order  when  n=0.  The 
equations  to  be  solved  in  this  study  will  be  of  sixth-order; 
that  is,  n=C.  Therefore,  six  solutions  will  be  computed 
numerically  and  linearly  combined  to  satisfy  the  six  boundary 
conditions.  All  of  the  derivations,  however,  will  be  done  for 
the  eighth-order  system  from  which  the  sixth-order  system  can 
easily  be  obtained  by  setting  n=0.  There  will  be  two  dif- 
ferent numerical  procedures  used:  one  will  be  used  for  the  neigh- 
borhood of  the  appex  where  difficulties  are  encountered  due  to 
small  angles,  and  another  for  the  rest  of  the  shell. 

4.4.1  Method  of  Integration  Away  From  the  Apex 

Included  in  the  total  solution  to  the  sixth-order 
system  of  shell  equations  are  four  non-growing  solutions  and 
one  growing  solution  in  each  direction.  One  solution  grows 
from  the  base  to  the  apex  and  the  other  grows  from  the  apex 
to  the  base.  These  two  growing  solutions  are  obtained  by 
numerically  integrating  the  sixth-order  shell  equations  from 
a point  near  the  apex  to  the  base  and  from  the  base  to  the 
same  point  near  the  apex.  The  four  non-growing  solutions  are 
determined  by  first  integrating  the  reduced  order  equations 
(B.3)  two  times  from  the  base  to  near  the  apex  and  two  times 
from  the  base  to  near  the  apex  and  two  times  from  near  the 
apex  to  the  base.  The  results  are  then  substituted  into 
Eqs . (B.l)  to  determine  the  displacements  and  stress  resul- 


tants. 
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4.4.2  Method  of  Integration  Near  the  Apex 

As  stated  earlier  the  singularity  at  the  apex,  due 
to  the  coordinate  system  chosen,  warrants  special  attention. 
The  problem  occurs  when  numerical  integration  is  carried 
out  in  the  neighborhood  of  the  apex.  Most  of  the 

coefficients  of  the  dependent  variables  are  divided  by 

. . 2 
sin  0 or  sin  0,  such  as 


n 

a sin  0 


u 


and 


2n2D ( 1-v) 

~J7  ~2~7 

a k sin  0 


Therefore,  these  terms  are  very  large  near  the  apex  and  un- 
defined at  0=0.  To  eliminate  this  problem  the  equations  of 
the  shell  in  the  neighborhood  of  the  apex  are  re-written  in 
terms  of  a new  independent  variable  Z,  defined  as  follows: 

z = Ln  (sin  0) 


It  is  seen,  that  as  sin  0,  and  thus  0,  approaches  zero,  z 
approaches  minus  infinity.  This  means  that  for  the  start  of 
integration  at  the  apex,  some  finite,  negative  number  for  z 
can  be  chosen  such  that  <|>  is  approximately  zero.  For  ex- 
ample, if  z=-10  then  0=4.5  x 10  \:ad.  *0  (see  Table  1).  The 
derivative  of  z with  respect  to  0 is 


dz  a cos  0 
3$  sin  0 
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A typical  derivative  term  in  the  shell  equations,  say  ^ is 
transformed  as  follows  when  the  independent  variable  is 
changed : 

du  _ du  , dz  _ du  . cos  d> 
d<J>  dz  dz  sin  $ 

All  of  the  derivative  terms  in  the  shell  equations  will  be  of 
the  form 


f*  (z) 


cos  <}> 
sin  <p 


To  solve  for  these  terms,  each  first-order  shell  equation  is 
divided  by  cos  <p/sin  <J>.  This  eliminates  all  of  the  1/sin  $ 
coefficients  in  the  shell  equations  which  means  a term  such 

as  nu/sin  <J>  becomes  nu/cos  - It  is  obvious,  that  1/cos  <J> 

. 2 
will  cause  no  difficulty  at  or  near  <p=0.  The  1/sin  <p  co- 
efficients do  not  cause  any  difficulty  at  (J>=0  either  because 
the  dependent  variables  by  which  they  are  multiplied  are  zero 
at  <J>=0  and  do  not  grow  rapidly  in  the  neighborhood  of  the 
apex.  The  interval  size  needed  starts  out  small,  but  can  be 
rapidly  increased  as  the  integration  moves  away  from  the  apex. 
Integration  continues  from  z=-10  to  z=-1.75  which  is  <j>=10°. 
Here,  the  derivatives  with  respect  to  z of  the  displacements 
and  stress  resultants  will  be  used  to  calculate  the  derivatives 
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with  respect  to  4> . At  this  point,  the  integration  will  con- 
tinue with  the  set  of  differential  Eqs . (4.3)  where  the 
independent  variable  is  <)>. 

When  integrating  the  reduced  order  equations  from 
the  base  to  the  apex,  again  at  <|)  = 10 0 derivatives  will  be 
calculated  in  terms  of  z.  From  this  point  to  the  apex  the 
differential  equations  written  in  terms  of  z will  be  inte- 
grated to  find  the  displacements  and  stress  resultants. 

In  this  interval  the  solutions  do  not  grow  to  the 
point  where  they  lose  their  independence.  Therefore,  from 
4> =0  to  <J>  = 10°  it  is  not  necessary  to  use  the  reduced  order 
equations  to  find  the  independent  solutions. 

4.5  Determination  of  the  Natural  Frequency 

In  determining  the  natural  frequencies  of  the  shell 
it  is  necessary  to  solve  the  homogeneous  portion  of  the  shell 
equations.  A combination  of  the  Holzer  method  (Holzer,  1921) 
and  a method  developed  by  Robinson  and  Harris  (Robinson  and 
Harris,  1971)  will  be  used  to  determine  the  natural  frequency. 
The  Holzer  method  will  be  used  to  obtain  a first  approximation 
for  a frequency,  then  the  Harris  method  will  be  used  to  improve 
upon  that  approximation. 

4.5.1  Holzer  Method 

The  application  of  the  Holzer  method  to  the  shell  is 
carried  out  as  follows: 


(1)  An  initial  frequency  oj  is  chosen  along 
with  six  independent  sets  of  initial  values  for 
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the  system  of  six  first-order  equations.  Six 
independent  initial  value  solutions  to  the 
homogeneous  equations  are  determined  using  the 
reduction  of  order  procedure. 

(2)  The  six  solutions  are  combined  to  satisfy 
the  six  homogeneous  boundary  conditions.  This 
is  only  possible  if  the  determinate  of  the  co- 
efficients, A ,is  zero.  (The  coefficients  in 
this  case  are  the  values  of  the  displacement 
and  the  stress  resultants  at  the  boundaries.) 
Therefore,  the  object  is  to  get  a zero  value 
of  the  determinant. 
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(3)  Values  of  w are  chosen  to  obtain  different 
values  of  A.  If  A changes  sign  we  know  that  a 
solution  to  the  homogeneous  system  of  equations 
lies  between  the  two  corresponding  values  of  w . 

The  value  of  oj  can  be  determined  only  within 
certain  limits  with  the  Holzer  method.  When  w gets 
within  a certain  range  of  the  true  w , the  deter- 
minent  of  the  coefficients  A begins  to  oscillate 
between  positive  and  negative  values  because  of  the 
small  differences  encountered  when  finding  the 
value  of  the  determinant.  Therefore,  it  becomes 
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necessary  to  go  to  a more  accurate  method 
such  as  the  Robinson-Harris  method  as  des- 
cribed in  the  next  section. 

4.5.2  Robinson-Harris  Method 

The  method  developed  by  Robinson  and  Harris  is 
an  application  of  the  Newton- Raphson  technique  of  improving 
an  eigenvalue  and  the  corresponding  eigenvector  from  approxi- 
mations of  the  eigenvalue  and  eigenvector.  In  this  section, 
the  Robinson-Harris  method  will  be  used  to  improve  the  ap- 
proximation of  the  natural  frequency  of  a spherical  shell. 

The  problem  to  be  solved  can  be  described  by 
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Equation  (4.7)  is  interperted  as  a non-linear  equa- 
tion in  A^  and  X^,  the  non-linear  term  being  A^BX^. 
Equation  (4.7)  is  linearized  about  the  reference  configuration 
corresponding  to  the  superscript  (i)  which  yields  in  linear- 
ized for 


A6X(l)  - A(l)B6X(l)  - 6A(l)BX(l)  = R(l)  (4.8) 

where  6X^  and  6A^  are  small  unknown  incremental  changes 
of  X ^ 1 ^ and  A^1^.  The  residual  R^  is  determined  from  equa- 
tion (4.7).  Equation  (4.8)  has  an  extra  unknown,  6A ^ , there- 
fore an  additional  equation  is  needed  to  get  a solution.  The 
additional  equation  is  taken  as 


o 


6X(1)  XU) 


dx  = 0 


(4.9) 


This  guarantees  that  the  allowable  change  in  the  vector  of  eigen- 
functions is  orthogonal  to  the  latest  vector  of  eigenfunctions,  to 
further  explain  this  procedure  the  actual  shell  equations  will  be 
used. 

The  shell  equations  are  first  put  into  the  form  of  Eqs . 

(4.6) , that  is , 

A | u,  v, w, s , ,M^  , , VQ 

Ab  ^u,v,w,s,I 
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where  A and  B are  8x8  matrices  composed  of  differential 
and  other  linear  operators.  An  approximate  eigenvalue  and 
eigenfunction  are  then  substituted  into  Equation  (4.10)  to 
give : 


(i)  (i)  (i)  (i)  „(i)' 

,ST  ,N<|>T  ' M4>T  ' V4»T  ,V6T  / 


/ (i)  (i)  (i] 

< UT  »VT  »WT 

' (i)  (i)  (i ) N(i)  M(i)  v(i)  v(i)\ 

T ,VT  ,wT  /N^  »M^t  ,V^t  ,Vgt  > 


,(i)  D(i)  D(i)  =(i)  D(i)  D(i)  B(i) 


Rv  ' Rv  ' R'  ' R'  ' Rl  ' R 

■5  i ^ a i i'r  , **£  , ^*7  , K 


8 


(4.11) 


The  approximate  eigenvalue  is  the  square  of  the  latest  frequency 
determined  in  the  Holzer  method.  The  approximate  eignefunction 
is  composed  of  stress  resultants  and  displacements  corresponding 
to  the  latest  frequency.*  Next,  Equations  (4.11)  are  linerized 
about  the  reference  configuration  yielding  the  following  non- 
homogeneous,  first-order  differential  equations: 
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dM4> » + f A (6u(l)  ,6v(l)  , 6w(l)  ,6s  (l)  J )=  -R*lJ 

id*1  6 V 4)  <J)  e y 6 

<t>  (1)  + f7^6u(l)  ,5v(l)  ,6w(l)  ,6s(l)  ,6M<Jl)  , 

6V0(l)  ,6X  (l)  ) = -R^l) 


ad  <f> 
d 


6 sin 


6 ii*{ve  sin  *>  ( 


1)  + f8^6u(l)  ,6vU)  ,6w(l)  ,6sU)  ,6N(J1)  ,6M(J1)  , 

(lM  = -R(i)  (4.12) 


vr'>)-  -r8 


We  now  have  eight  equations  in  nine  unknowns,  the  unknowns  being 
the  incremental  changes  in  the  displacements  and  stress  resul- 
tants, and  the  incremental  change  in  the  square  of  the  natural 
frequency.  Again,  as  in  Equation  (4.8),  we  have  one  more  un- 
known than  equations;  thus,  an  additional  side  condition  is 
needed  to  solve  the  above  system  of  differential  equations.  In 
this  case,  the  side  condition  is  taken  as 


4 2 it  **  ’ 

f 0 f pha2  sin  <p  ( u(l)6u(l)  + w(l)<5w(l)  + v(l)6v(l)) 


d 9d <J>  =>  0 


which  simplifies  to 


♦ 

- 

* 

1 

f .l« 

u(ll«u(i)  ♦ w(i,«wa)  + v(i)Sv(i) 

1 

d<p  = 0 

(4.13) 
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This  is  the  equivalent  to  saying  that  the  incremental  change 
in  the  displacement  vector  is  orthogonal  to  the  latest  dis- 
placement vector  with  respect  to  the  mass.  We  now  have  nine 
equations  in  nine  unknowns. 

The  boundary  value  problem  described  by  Equations 
(4.12)  and  (4.13)  is  solved  by  the  method  of  reduction  of 
order  of  a system  of  differential  equations  as  described  in 
Chapter  2,  with  the  exception  that  for  Equation  (4.13)  there 
is  no  reduction  of  order  just  quadrature  by  Simpson's  Rule. 
The  initial  values  used  in  the  solution  of  the  boundary 
value  problem  are  listed  in  Table  2 . Ten  independent  solu- 
tions, nine  homogeneous  and  one  particular  are  determined 
using  the  initial  values  in  the  reduction  of  order  method. 

The  solutions  are  then  combined  to  satisfy  the 
following  eight  homogeneous  boundary  conditions*: 


0 = 6u1Q(0)  + j^XjfiujiO) 


(u=0  at  <|>=0) 


0 = 6wi0(0)  + ^1-,x46w^(0) 


j=l  j j 


(w=0  at  <|>=0) 


0 =6v1q(0)  + j|1Xj6vj(0) 


(v=0  at  <f>=0) 


0 = 6s1q(0)  + jliXjdSjto) 


(s=0  at  <p  =0 ) 


For  n=0  or  n=l  the  boundary  conditions  would  be  changed 
appropriately . 
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(w=0  at  <p=<pQ) 


(v=0  at  <|>=4>o) 

(s*°  at  ♦‘♦o’  (4.14) 


where  j is  the  initial  value  solution  number  and  x^  is  the 
weighting  coefficient.  The  additional  equation  needed  to  deter- 
mine the  weighting  coefficients  is  one  involving  the  side  condi- 
tion and  expressed  as 


0 = f si 

r\  ~ 


in  cj) 


u^6u,„  + w^6w,„  + v^6v 


*10 


'10 


10 


d<f> 


(7  V 

+ . Z.  x . f 

3=1  3 I 

0 J 


0 

sin  4> 


w(i^6w..  + v(l)6Vj 
(4.15) 


d<J> 


From  Equations  (4.14)  and  (4.15)  the  x^'s  are  calculated. 
Once  their  values  are  known  the  total  incremental  change  in 
eigenvalue  and  displacements  and  stress  resultants  can  be  com- 
puted, that  is. 


6u 


(k) 


= 6u 


(k) 

10 


j=lXj 


6u 


(k) 

j 


The  other  variables  and  their  derivatives  can  be  similarly  ex- 
pressed. The  incremental  changes  are  then  added  to  the  total 
solution  to  get 

(k+1)  (k)  . , (k) 

u — u + 6u 

The  new  reference  position  is,  therefore,  the  (k+l)th  ap- 
proximation. This  n^w  approximation  is  then  substituted  back 
into  Equation  (4.11)  to  compute  new  residuals.  The  new  resi- 
duals are  checked  to  see  if  they  are  within  a certain  toler- 
ance. If  they  are,  the  procedure  is  terminated.  If  not,  the 
complete  process  is  repeated  and  new  residuals  are  calculated 
and  checked  for  accuracy. 
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5.  NUMERICAL  RESULTS  OF  THE  APPLICATION  OF  THE  METHOD  OF 

REDUCTION  OF  ORDER 


5. 1 General  Remarks 

In  this  Chapter,  the  results  of  the  application  of 
the  proposed  method  to  a system  with  multiple  edge  effects  is 
presented.  In  addition,  the  results  of  the  calculations  of 
free  vibration  of  spherical  shell  segments  is  also  present- 
ed. The  multiple  edge  effect  example  is  solved  using  the 
methods  presented  in  Chapters  2 and  3,  and  the  shell  example 
is  solved  using  the  methods  presented  in  Chapters  2 and  4 . 

The  main  objective  of  presenting  these  examples  is  to  demon- 
strate some  of  the  possible  uses  of  the  proposed  method. 

The  multiple  edge  effect  example, 

wlv  - (kj+k^w"  + k^k^w  = p ( x ) (3.2) 

given  in  Section  5.2  contains  two  rapidly  growing  solutions 
which  are  effectively  separated  numerically  by  the  proposed 
method.  The  results  of  the  calculations  are  compared  with 
the  exact  solution  of  Equation  (3.2). 

The  examples  given  in  Section  5.3  are  spherical 
shell  segments  where  the  natural  frequencies  and  mode  shapes 
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have  been  computed  for  <f>Q  = 5°,  10°,  18°.  The  natural  fre- 
quencies and  the  mode  shapes  are  found  for  these  shell  seg- 
ments using  the  Holzer  method  and  the  method  developed  by 
Robinson  and  Harris  as  presented  in  Chapter  4.  In  the  appli- 
cation of  the  Holzer,  and  Robinson  and  Harris  methods,  it  is 
necessary  to  repeatly  solve  two-point  boundary  value  problems . 
The  reduction  of  order  method  is  used  to  solve  these  boundary 
value  problems. 

5 . 2 Example  With  Multiple  Edge  Effects 

In  solving  a problem  with  multiple  edge  effects,  the 
difficulty  lies  in  first  trying  to  seperate  the  two  growing 
solutions  and  then  trying  to  obtain  a non-growing  particular 
solution.  The  example  described  by  Equation  (3.2)  was  first 
re-written  as  four  first  order  equations,  as  described  in 
Chapter  3,  for  conveninece.  The  equations  are 

t'  - (k^  + k^) s + k^k^q  = p(x) 

q'  - r = 0 

r'  - s = 0 


s'  - t = 0 


--  - - 


*2  = 7 


The  interval  length  was  1 and  p(x)  was  10,000.  Since  the 

lex 

solutions  were  of  the  form  e , it  is  obvious,  that  both 

1 ft  *7  V 

solutions,  w1=e  and  w2=e  , grow  rapidly.  The  boundary 

conditions  used  are  as  follows: 


The  values  of  k were 


k1  = 10 


w (0 ; =1  w'  (L)  = -4 


w"  (L)  = 0 w"'  (0  ) = 0 


(5.1) 


The  exact  solution  was  determined  using  Equation  (3.1) 
If  this  equation  is  differentiated  three  times  we  have 


k.x  k2x  -k.x  -k2x 

w = C1e  1 + C2e  + C3e  + C4e 


k.x  k_x  -k.x  -k2x 

w'  = k1C1e  1 + k2C2e  z - k-^e  1 - k2C4e 


o k.x  k x 2 —k.x  2 — k.x 

w"  = k^Cxe  1 + k2C2e  1 + k^e  x + k2C4e 


, k.x  k_x  , k.x  , -k_x 

w'"  = k^e  + k2C2e  d - kjc3e  1 - k^e 


(5.2) 
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This  is  the  homogeneous  solution  of  Eq.  (3.1).  The  particular 
solution  is 


wp  = p (x) / (k^+k^) 

The  coefficients  in  Eq.  (5.2)  are  determined  by  setting  the 
homogeneous  and  particular  solutions  and  their  derivatives 
equal  to  boundary  conditions  (5.1).  Thus,  once  the  coefficients 
are  determined,  the  exact  solution  can  be  calculated. 

5 . 3 Numerical  Results  for  Multiple  Edge  Effects  Example 

Data  given  in  Table  3 shows  the  rate  of  growth  of  the 
fastest  growing  solution.  This  solution  was  obtained  by  straight- 
forward numerical  integration  of  the  homogeneous  part  of  Eq. 

(3.2).  It  is  seen  that  when  x is  at  all  large 

10^  = 102w^  = 10w£  = w" 

This  is  exactly  what  we  have  in  the  exact  homogeneous  solution, 
that  is. 
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w = e 


k^x 


w'  = k.e 


k^x 


w”  = k7e 


2 kl* 


3 klx 

w'"=  kje  l 


The  above  is  also  true  for  the  solution,  w^,  which  grows 
from  L to  0. 

Data  given  in  Table  4 is  the  slower  growing  solu- 
tions. This  solution  was  obtained  using  Eqs . (3.6)  and  Eqs . 
(3.7).  It  is  seen  that 


w„  = e 


k2x 


w. 


k2e 


k2X 


w. 


„2  k X 
k2e 


Here  we  have 


w. 


3 k2X 

k2e 


where  k2  = 7 


73 * * * 7w2  = 72w^  = 7w"  = w"’ 
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The  above  is  also  true  for  the  slower  growing  solution,  w^  , 
which  grows  from  L to  0 . 

From  experience  it  has  been  found  that  a value  of 
c of  about  three  times  the  maximum  value  of  w^  yields  the 
best  results.  The  value  of  c should  be  chosen  such  that 
the  quantity  q^+c  does  not  vary  greatly.  A value  of  c 
much  larger  than  just  what  is  needed  to  make  (q^+c)  con- 
stant increases  round-off  error  and  decreases  the  accuracy 
of  the  results. 

The  data  given  in  Tables  3 and  4 clearly  indicate 
that  the  two  growing  solutions  from  each  end  of  the  interval 
are  accurately  separated.  Thus,  if  they  are  accurately 
separated  they  then  can  be  linearly  combined  when  determining 
the  total  solution  to  Eq.  (3.1). 

The  non-growing  particular  solutions  listed  in 
Tables  5 and  6 were  determined  using  Eqs.  (3.8)  and 
(3.10),  respectively.  Both  tables  list  particular  solutions 
to  Eq . (3.6). 

In  Table  7 total  solution  number  1,  containing  parti- 
cular solution  1 and  total  solution  number  2,  containing  parti- 
cular solution  2 are  listed.  As  can  be  seen,  when  compared 
to  the  exact  solution,  total  solution  number  1 is  more  ac- 
curate when  x is  small  and  total  solution  2 is  more  accurate 
when  x is  large.  Also,  it  is  seen  that  solutions  1 and  2 are 
almost  identical  in  the  middle  of  the  interval  and  both  compare 
very  well  with  the  exact  solution  at  this  point. 
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Table  8 shows  comparisons  between  the  exact  so  - 
lution  and  the  solution  obtained  by  the  reduction  of  order 
method.  The  solutions  obtained  by  reduction  of  order  con- 
tain portions  of  two  different  solutions  as  described 
above . 

The  calculations  for  the  reduction  of  order  method 
were  done  on  an  IBM  360  computer,  using  single  precision. 

The  exact  solution  was  calculated  on  the  same  computer,  but 
using  double  precision. 

5 . 4 Natural  Frequencies  and  Mode  Shapes  of  Spherical  Shell 
Segments 

The  numerical  results  of  the  method  of  analysis  used 
in  the  determination  of  natural  frequencies  and  mode  shapes  of 
spherical  shells  is  presented  in  this  section.  The  Holzer 
method,  and  the  Robinson-Harris  method  are  used  to  deter- 
mine the  natural  frequencies  and  mode  shapes  of  the  spherical 
shell  segments.  The  reduction  of  order  method  is  used  to  solve 
the  two-point  boundary  value  problems  which  occur  in  both 
methods. 

The  shell  segments  used  in  this  study  had  opening 
angles  of  5°,  10°,  and  18°.  The  h/a  ratio  was. 1/100  and  the 
base  of  the  shell  was  clamped,  that  is, 

u = w = s = 0 (s  = dw/d<(>) 
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The  Holzer  method  is  used  to  determine  an  approxi- 
mate value  of  the  natural  frequency  parameter 

w = < V-  yE"T 

wo  ° V pa2  (i-v2) 

This  approximation  is  improved  upon  by  the  Robinson-Harris 
method.  Table  9 gives  a comparison  between  the  natural  fre- 
quency parameter  o determined  by  the  two  methods. 

When  two  successive  values  of  w are  within  a cer- 
tain tolerance,  the  last  value  of  w determined  in  the  Hol- 

_2 

zer  method  becomes  the  first  approximation  of  X(X=w  ) in 
the  Robinson-Harris  method.  Table  10  lists  the  values  of  u 
determined  by  the  Robinson-Harris  method  from  the  first  ap- 
proximation to  the  final  value  and  the  values  of  the  resi- 
duals R^.  As  can  he  seen  from  the  table,  there  is  a very 
rapid  convergence  when  starting  with  a fairly  good  approxi- 
mation of  X. 

The  problem  at  the  apex,  due  to  the  coordinate  system 
chosen,  is  treated  by  changing  the  indpendent  variable  from 
<P  to  z as  described  in  Section  4.4.2.  That  is,  z=ln 
(sin  4>)  . By  changing  the  independent  variable  to  z the 

_ 5 

integration  can  begin  very  close  to  <p=0  (<j>=4.5  x 10  rad.) 
without  difficulty.  The  size  of  the  intervals  is  very  small 
at  4>=0,  but  is  increased  rapidly  as  4>  increases.  Table  11 
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shows  how  the  intervals  are  increased  as  4>  increases  and 
also  shows  the  values  of  the  displacements. 

A comparison  is  made  between  the  method  applied  in 
this  study  and  work  done  by  Zarghamee  in  Table  12.  It  can  be 
observed  that  there  is  a fairly  good  comparison  between  the 
values  of  the  natural  frequency  parameter  u) . 

11 
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6.  GENERAL  CONCLUSIONS 

6 . 1 Multiple  Edge  Effects 

The  reduction  of  order  method  developed  in  this 
study  provides  a procedure  for  determining  accurate  solutions 
of  two-point  boundary  value  problems  when  edge  effects  are 
present.  Extremely  good  agreement  was  obtained  when  the 
exact  solution  was  compared  to  the  reduction  of  order  solu- 
tion of  the  multiple  edge  effect  example.  Although  this 
method  was  only  applied  to  examples  of  problems  with  one  and 
two  edge  effects  it  can  be  extended  to  problems  where  there 
are  more  than  two  edge  effects.  A general  procedure  was 
described  for  the  solution  of  an  nth-order  system  with  pos- 
sibly n/2  growing  solutions.  The  procedure  works  well  for 
homogeneous  as  well  as  non-homogeneous  problems. 

6 . 2 Free  Vibrations  of  Spherical  Shells 

The  method  of  analysis  used  in  this  study  for  deter- 
mining the  natural  frequencies  and  mode  shapes  of  spherical 
shells  was  a combination  of  the  Holzer  and  Robinson-Harris 
method.  The  reduction  of  order  method  was  used  to  solve  the 
two-point  boundary  value  problems  which  occurred  in  these 
methods.  The  difficulties  created  at  the  apex  by  the  coor- 
dinate system  were  effectively  dealt  with  by  changing  the  in- 
dependent variable  in  the  shell  differential  equations.  This 
allowed  numerical  integration  to  begin  at  a point  close  enough 
to  4)  = 0 such  that  for  all  practical  purposes  <J>  could  be 
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TABLE  1.  VALUES  OF  THE  QUANTITY  z, 
THE  COLATITUDE  <J>  AND  SIN  (J> 

Z sin  $ (f>  (rad.) 


— 00 

0 

0 

-10 

4.5 

X 

10-5 

4.5 

X 

10~5 

-5 

6.7 

X 

10~3 

6.7 

X 

10"3 

-4 

1.8 

X 

ID"2 

00 

• 

H 

X 

(N 

1 

O 

H 

-3 

4.98 

X 

IQ"2 

4.98 

X 

CM 

1 

O 

H 

“2  .1353  .1353 
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TABLE  3.  FASTEST  GROWING  SOLUTION  FROM  0 TO  L 


w 


1 


w: 


w: 


! 

i 

I 


1 

0.0 

0.0 

2 

0 . 4582567E-12 

0 . 3078551E-10 

3 

0 . 21224 14E-11 

0 . 7007898E-10 

4 

0. 5274820E-11 

0.1209362E-09 

5 

0 . 1042904E-10 

0 . 1912656E-09 

6 

0 . 1840207E-10 

0 . 2914942E-09 

7 

0 . 304 1961E- 10 

0.4358829E-09 

8 

0.4827243E-10 

0.6442988E-09 

9 

0 . 74  542  30E-10 

0 . 944654  IE-09 

10 

0 . 1129249E-09 

0 . 1376300E-08 

11 

0 . 1686872E-09 

0.1994767E-08 

12 

0 • 2493132E-09 

0 . 2878414E-08 

13 

0.3654148E-09 

0 . 4137735E-08 

14 

0 . 5320133E-09 

0.5928442E-08 

15 

0.7703413E-09 

0.8469794E-08 

16 

0 . 1110374E-08 

0.1207026E-07 

17 

0 . 1594377E-08 

0 . 1716359E-97 

18 

0. 2281902E-08 

0.2435923E-07 

19 

0 . 3256774E-08 

0.3451311E-07 

20 

0 . 4636888E-08 

0. 4882660E-07 

21 

0 . 6587982E-08 

0. 6898506E-07 

22 

0 . 9342873E-08 

0 . 973524  3E-07 

23 

0 . 1322844E-07 

0 . 1372426E-06 

24 

0 . 1870341E-07 

0.1933000E-06 

25 

0 . 2641126E-07 

0.2720309E-06 

26 

0 . 3725426E-07 

0 . 3825495E-06 

27 

0 . 5249710E-07 

0. 5376176E-06 

28 

0. 7391185E-07 

0 . 7551022E-06 

29 

0 . 103  9811E-0  6 

0 . 1060011E-05 

30 

0 . 1461816E-06 

0.1487343E-05 

31 

0 . 2053816E-06 

0 . 2086076E-05 

0.9999999E-09 
0 . 1068107E-08 
0.1312822E-08 
0 • 1766964E-08 
0 . 2489540E-08 
0 . 3573659E-08 
0 . 5157670E-08 
0 . 7441709E-08 
0 . 1071142E-07 
0.1537175E-07 
0 . 2199467E-07 
0.3138624E-07 
0 . 4468085E-07 
0. 6347335E-07 
0 . 9000433E-07 
0.1274198E-06 
0 . 1801344E-06 
0 . 2543407E-06 
0.3587215E-06 
0 . 5054470E-06 
0 . 7115705E-06 
0.1000979E-05 
0.1407127E-05 
0 . I976853E-05 
0 . 2775732E-05 
0.3895539E-05 
0 . 5464695E-05 
0 . 7662891E-05 
0 . 1074151E-04 
0 . 1505218E-04 
0 . 2108664E-04 


0.0 

0 . 4564470E-08 
0.1022545E-07 
0.1721998E-07 
0 . 2644401E-07 
0.3906096E-07 
0 . 5663972E-07 
0.8133958E-07 
0.1161645E-06 
0.1653132E-06 
0 . 2346686E-06 
0.3324797E-06 
0 . 4703170E-06 
0.6644074E-06 
0 . 9375022E-06 
0.1321489E-05 
0 . 1861063E-05 
0 . 2618828E-05 
0 . 3682509E-05 
0 . 5174872E-05 
0 . 7267845E-05 
0 . 1020206E-04 
0 . 1431426E-04 
0.2007559E-04 
0.2814534E-04 
0.3944570E-04 
0 . 5526657E-04 
0.7741200E-04 
0.1084048E-03 
0 . 1517726E-03 
0 . 2124478E-03 
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TABLE  4 . SLOWER  GROWING  SOLUTION  FROM  0 TO  L 


w2  W2  W2  W2 

0 0.0  0 . 9999999E-09  0 . 4893369E-08  -0 . 9999999E-09 

1 0.3233699E-10  0 . 1150013E-08  0 . 5173288E-08  0 . 1902207E-07 

2 0. 7334304E-10  0 . 1337381E-08  0 . 6137736E-08  0 . 3832426E-07 

3 0 . 1222734E-09  0 . 1567526E-08  0 . /725976E-08  0 . 5665011E-07 

4 0 . 1794046E-09  0 . 1861117E-08  0 . 9936436E-08  0 . 7580849E-07 

5 0 . 2477141E-09  0 . 2239733E-08  0 . 1282402E-07  0 . 97396I4E-07 

6 0. 3304412E-09  0 . 2727664E-08  0. 1649552E-07  0 . 1229456E-06 

7 0 . 4 317087E-09  0 . 3353642E-08  0 . 2110981E-07  0 . 1540483E-06 

8 0.5566982E-09  0 . 4152621E-08  0 . 2688151E-07  0 . 1924676E-06 

9 0.7118923E-09  0 . 5167777E-08  0 . 3408834E-07  0 . 2402354E-06 

10  0 . 9053860E-09  0 . 6452762E-08  0 . 4308220E-07  0 . 2997586E-06 

11  0 . 1147280E-08  0 . 8074451E-08  0 . 54 30320E-07  0 . 3739297E-06 

12  0 . 1450176E-08  0 . 1011610E-07  0 . 6829748E-07  0 . 4662297E-06 

13  0 . 1829791E-08  0 . 1268134E-07  0 . 8573915E-07  0 . 5808624E-06 

14  0 . 2305712E-08  0 . 1589892E-07  0 . 1074567E-06  0 . 7228800E-06 

15  0 . 2902335E-08  0 . 1992818E-07  0 . 133464 0E-06  0 . 8983016E-06 

16  0.3649991E-08  0 . 2496602E-07  0 . 1679959E-06  0 . 1114311E-05 

17  0.4586326E-08  0 . 3125524E-07  0 . 2095467E-06  0 . 1379308E-05 

18  0 . 5757990E-08  0 . 3909298E-07  0 . 2609177E-06  0 . 1702996E-05 

19  0 . 7222596E-08  0 . 4884231E-07  0 . 3242503E-06  0 . 2096480E-05 

20  0 . 9051252E-08  0 . 6094376E-07  0 . 4 020740E-06  0 . 2571745E-05 

21  0 . 1133129E-07  0 . 7593042E-07  0 . 4973326E-06  0 . 3141613E-05 

22  0 . 14 1696bE-07  0 . 9444091E-07  0 . 6134014E-06  0 . 3818728E-05 

23  0 . 1769661E-07  0. 1172341E-06  0 . 7540557E-06  0 . 4614005E-05 

24  0 . 2206974E-07  0 . 1452017E-06  0 . 9234027E-06  0 . 5534937E-05 

25  0 . 274 7902E-07  0 . 1793760E-06  0 . 1125707E-05  0 . 6581395E-05 

26  0 . 3415118E-07  0 . 2209254E-06  0 . 1364999E-05  0 . 7740222E-05 

27  0. 4235375E-07  0. 27114 31E-06  0 . 1644526E-05  0 . 8977018E-05 

28  0. 5239974E-07  0 . 3314017E-06  0 . 1965818E-05  0 . 1022033E-04 

29  0 . 6464904E-07  0 . 4030799E-06  0 . 2 327171E-05  0 . 1134491E-04 

30  0 . 7950808E-07  0 . 4874310E-06  0 . 2721223E-05  0 . 1214235E-04 
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TABLE  5.  PARTICULAR  SOLUTION,  W , 

pi 


w 

w' 

w" 

w"' 

P2 

P2 

P3 

P4 

0 

0.0 

0.1097E 

02 

-0.4511E 

02 

-0.1097E 

02 

1 

0.3099E 

00 

0.9624E 

01 

-0.4419E 

02 

0.6559E 

02 

2 

0.6076E 

00 

0.8193E 

01 

-U.4114E 

02 

0.1110E 

03 

3 

0.8593E 

00 

0.6884E 

01 

-0.3706E 

02 

0.1295E 

03 

4 

0.1070E 

01 

0.5718E 

01 

-0.3266E 

02 

0.1317E 

03 

5 

0.1243E 

01 

0.4698E 

01 

-0.2836E 

02 

0.1242E 

03 

6 

0.1385E 

01 

0.3817E 

01 

-0.2441E 

02 

0.1115E 

03 

7 

0.1500E 

01 

0.3060E 

01 

-0.2093E 

02 

0.9645E 

02 

8 

0.1591E 

01 

0.2411E 

01 

-0.1796E 

02 

0.8059E 

02 

9 

0.1662E 

01 

0.1853E 

01 

-0.1553E 

02 

0.6492E 

02 

10 

0.1716E 

01 

0.1366E 

01 

-0.1362E 

02 

0.4988E 

02 

11 

0.1754E 

01 

0.9365E 

00 

-0.1218E 

02 

0.3571E 

02 

12 

0.1779E 

01 

0.5464E 

00 

-0.1121E 

02 

0.2240E 

02 

13 

0.1791E 

01 

0.1810E 

00 

-0.1068E 

02 

0. 9793E 

01 

14 

0.1791E 

01 

0.1738E 

00 

-0.1056E 

02 

-9.2391E 

01 

15 

0.1780E 

01 

-0.5305E 

00 

-0.1083E 

02 

-0.1379E 

02 

16 

0.1756E 

01 

-0.9031E 

00 

-0.1147E 

02 

-0.2472E 

02 

17 

0.1719E 

01 

-0.1298E 

01 

-0.1248E 

02 

-0.3482E 

02 

18 

0.1668E 

01 

-0.1733E 

01 

-0.1375E 

02 

-0.4348E 

02 

19 

0.1602E 

01 

-0.2218E 

01 

-C.1522E 

02 

-0.4823E 

02 

20 

0.1519E 

01 

-0.2759E 

01 

-0.1685E 

02 

-0.4879E 

02 

21 

0.1417E 

01 

-0.3357E 

01 

-0.1856E 

02 

-0.4321E 

02 

22 

0.1295E 

01 

-0.4004E 

01 

-0.1997E 

02 

-0.3107E 

02 

23 

0.1153E 

01 

-0.4667E 

01 

-0.2046E 

02 

0.1367E 

00 

24 

0.9852E 

00 

-0. 5336E 

01 

-0.1977E 

02 

0.5426E 

02 

25 

0.7903E 

00 

-0.5889E 

01 

-0.1610E 

02 

0.1528E 

03 

26 

0.5767E 

00 

-0.6255E 

01 

-0.7973E 

01 

0.3146E 

03 

27 

0.3672E 

00 

-0.6219E 

01 

0.8320E 

01 

0.5879E 

03 

28 

0.1724E 

00 

-0.5641E 

01 

0.3617E 

02 

0.9978E 

03 

29 

0.2710E- 

■01 

-0.3688E 

01 

0.7590E 

02 

0 . 161 IE 

04 

30 

0.9033E- 

-02 

0.2148E 

00 

0.1461E 

03 

0.2464E 

04 

98 


0 

0.9033E- 

-02 

-0.2148E 

00 

0.1461E 

03 

-0.2464E 

04 

1 

0.1709E 

01 

0.3457E 

01 

0.8196E 

02 

-0.1684E 

04 

2 

0.1565E 

00 

0.5508E 

01 

0.3883E 

02 

-0.1047E 

04 

3 

0.3452E 

00 

0.6297E 

01 

0.9965E 

01 

-0.6186E 

03 

4 

0.5544E 

00 

0.6267E 

01 

0.6617E 

01 

-0.3392E 

03 

5 

0.7693E 

00 

0.5944E 

01 

-0.1548E 

02 

-0.1655E 

03 

6 

0 . 967 IE 

00 

0.5396E 

01 

-0.1955E 

02 

-0.6207E 

02 

7 

0.1137E 

01 

0.4747E 

01 

-0.2044E 

02 

-0.4578E 

01 

8 

0. 1282E 

01 

0.4072E 

01 

-0.1999E 

02 

0.2909E 

02 

9 

0.1406E 

01 

0.3419E 

01 

-0.1869E 

02 

0.4249E 

02 

10 

0.1510E 

01 

0.2818E 

01 

-0.1706E 

02 

0.4863E 

02 

11 

0.1594E 

01 

0.2269E 

01 

-0.1536E 

02 

0.4848E 

02 

12 

0.1662E 

01 

0.1789E 

01 

-0.1390E 

02 

0.4397E 

02 

13 

0.1715E 

01 

0.1341E 

01 

-0.1260E 

02 

0.3587E 

02 

14 

0.1753E 

01 

0.9407E 

00 

-0.1156E 

02 

0.2575E 

02 

15 

0.1778E 

01 

0.5675E 

00 

-0.1087E 

02 

0.1495E 

02 

16 

0.1791E 

01 

0.2087E 

00 

-0.1056E 

02 

0.3559E 

01 

17 

0.1792E 

01 

-0.1455E 

00 

-0.1064E 

02 

-0.8637E 

01 

18 

0.1781E 

01 

-0.5089E 

00 

-0.1114E 

02 

-0.2111E 

02 

19 

0.17575 

01 

“0 . 896  IE 

00 

-0.1207E 

02 

-0.3435E 

02 

20 

0 . 17  2 IE 

01 

-0.1321E 

01 

-0.1345E 

02 

-0.4841E 

02 

21 

0.1669E 

01 

-0.1801E 

01 

-0.1532E 

02 

-0.6338E 

02 

22 

0.1599E 

01 

-0.2352E 

01 

-9.1770E 

02 

-0.7901E 

02 

23 

0.1510E 

01 

-0.2991E 

01 

-0.2061E 

02 

-0.9487E 

02 

24 

0.1398E 

01 

-0. 3736E 

01 

-0.2404E 

02 

-0.1101E 

03 

25 

0.1259E 

01 

-0.4604E 

01 

-0.2795E 

02 

-0 . 12  3 IE 

03 

26 

0.1089E 

01 

-0.5609E 

01 

-0.3222E 

02 

-0.1313E 

03 

27 

0.8822E 

00 

-0.6761E 

01 

-0.3663E 

02 

-0.1304E 

03 

28 

0.6348E 

00 

-0.8056E 

01 

-0.4076E 

02 

-0.1138E 

03 

29 

0.3419E 

00 

-0.9477E 

01 

-0.4396E 

02 

-0.7173E 

02 

30 

0.0 

-0.1097E 

02 

-0.4511E 

02 

0.1097E 

02 
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TABLE  7.  COMPARISON  BETWEEN  TOTAL  SOLUTIONS 
1 AND  2 AND  THE  EXACT  SOLUTION 


w_ 

w_, 

EXACT 

T1 

T2 

SOLUTION 

0 

0.9999999E 

00 

0 . 1000000E 

01 

0 . 1000000D 

01 

1 

0 . 1159199E 

00 

0.11185422 

01 

0 . 1159139D 

01 

2 

0 . 1311897E 

01 

0 . 1267189E 

01 

0 . 1311907D 

01 

3 

0 . 1441156E 

01 

0.1395608E 

01 

0 . 1441141D 

01 

4 

0 . 1549276E 

01 

0.150880 IE 

01 

0 . 1549242D 

01 

5 

0 . 1638834E 

01 

0 . 1613165E 

01 

0 . 1638786D 

01 

6 

0 . 171234  3E 

01 

0 . 1698983E 

01 

0.1712292D 

01 

7 

0 . 1772130E 

01 

0 . 1762797E 

01 

0 . 177208  3D 

01 

8 

0 . 1820264E 

01 

0 . 1811711E 

01 

0 . 1820218D 

01 

9 

0.1858508E 

01 

0 . 1851007E 

01 

0.1858473D 

01 

10 

0.1888365E 

01 

0 . 1882613E 

01 

0 . 1888334D 

01 

11 

0.1910976E 

01 

0. 1906458E 

01 

0 . 1911009D 

01 

12 

0.1927374E 

01 

0 . 1923646E 

01 

0 . 1927449D 

01 

13 

0.1938210E 

01 

0.1935847E 

01 

0 . 1938364D 

01 

14 

0.1943997E 

01 

0.1942519E 

01 

0 . 1944242D 

01 

15 

0 . 1944908E 

01 

0 . 1943964E 

01 

0 . 1945359D 

01 

16 

0 . 1941240E 

01 

0.1940779E 

01 

0 . 1941795D 

01 

17 

0 . 19  32470E 

01 

0 . 19  3268  5E 

01 

0 . 1933434D 

01 

18 

0.1918213E 

01 

0 . 1919416E 

01 

0 . 1919973D 

01 

19 

0.1898513E 

01 

0 . 1900434E 

01 

0 . 1900917D 

01 

20 

0.1872685E 

01 

0.1875299E 

01 

0 . 1875584D 

01 

21 

0.1839049E 

01 

0 . 1842897E 

01 

0.1843108D 

01 

22 

0.1798345E 

01 

0 . 1802300E 

01 

0 . 1802448D 

01 

23 

0. 1748493E 

01 

0 . 1752319E 

01 

0 . 1752421D 

01 

24 

0.1684845E 

01 

0.1691689E 

01 

0 . 1691750D 

01 

25 

0 . 1602128E 

01 

0 . 1619122E 

01 

0.1619162D 

01 

26 

0.1503093E 

01 

0 . 1533532E 

01 

0 . 1533556D 

01 

27 

0 . 1400103E 

01 

0 . 14  34243E 

01 

0 . 1434261D 

01 

28 

0.1286058E 

01 

0.1321447E 

01 

0 . 1321465D 

ei 

29 

0. 1167562E 

01 

0 . 1196843E 

01 

0 . 1196354D 

01 

30 

0 . 1078046E 

01 

0 . 1064647E 

01 

0 . 1064606D 

01 

100 


TABLE  8.  SOLUTION  OBTAINED  BY  REDUCTION  OF 
ORDER  COMPARED  WITH  THE  EXACT  SOLUTION 

REDUCTION 

EXACT  OF 

SOLUTION  ORDER 


0 

0 . 100000D 

01 

0 . 100000E 

01 

1 

0.115914D 

01 

0 . 1159  HE 

01 

2 

0 . 131191D 

01 

0 . 131190E 

01 

3 

0 . 144114D 

01 

0 . 144116E 

01 

4 

0 . 154924D 

01 

0.154928E 

01 

5 

0.163879D 

01 

0 . 163883E 

01 

6 

0 . 171229D 

01 

0 . 171234E 

01 

7 

0 . 177208D 

01 

0 . 177213E 

01 

8 

0.182022D 

01 

0 . 182026E 

01 

9 
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01 

0 . 185851E 

01 

10 

0 . 188833D 

01 

0 . 188837E 

01 

11 
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01 

0.191098E 

01 

12 
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01 

0 . 192737E 

01 

13 
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01 
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01 

14 
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01 
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01 

15 
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01 

0 . 194491E 

01 

16 
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01 
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01 

17 

0.193343D 

01 

0.193269E 

01 

18 

0 . 191997D 

01 

0.191942E 

01 

19 

0 . 190092D 

01 

0.190048E 

01 

20 

0 . 187558D 

01 

0 . 187530E 

01 

21 

0.184310D 

01 

0 . 184290E 

01 

22 

0 . 180245D 

01 

0.180230E 

01 

23 

0 • 175242D 

01 

0. 175232E 

01 

24 

0 . 169175D 

01 

0 . 169169E 

01 

25 

0 . 161916D 

01 

0.161912E 

01 

26 

0 . 153356D 

01 

0 . 153353E 

01 

27 

0 . 143426D 

01 

0.143424E 

01 

28 

0 . 132146D 

01 

0 . 132145E 

01 

29 

0 . 119685D 

01 

0 . 119684E 

01 

30 

0 . 106461D 

01 

0 . 106465E 

01 
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TABLE  12.  NATURAL  FREQUENCY  PARAMETER,  u> 
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FIG.  6.  MOMENT  STRESS  RESULTANTS  ACTING  ON  TWO 
ADJACENT  ELEMENTS  OF  A SHELL  SEGMENT. 
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FIG.  7 . COORDINATE  SYSTEM  USED  FOR  DERIVATION  OF 
KIRCHHOFF  FORCE . 


'’A* 
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FIG.  8.  EQUIVALENT  FORCE  SYSTEM  ACTING  ON  TWO 
ADJACENT  ELEMENTS  OF  A SHELL  SEGMENT 
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APPENDIX  A 

DERIVATION  OF  A GENERAL  SYSTEM  OF  FIRST-ORDER 
EQUATIONS  OF  MOTION  FOR  A SPHERICAL  SHELL 

A. 1 Introduction 

In  this  appendix,  eight  first-order  differential 
equations  of  motion  for  an  elastic  spherical  shell  will  be 
derived.  This  will  be  accomplished  by  considering  the 
equations  of  motion,  the  strain-displacement  relationships 
and  the  relationship  between  stress  and  strain  of  a typical 
shell  element. 

Throughout  the  derivation,  Love's  first  approxima- 
tion (Love,  1881)  will  apply.  That  is,  deformations  shall  be 
small,  a condition  of  plain  stress  shall  exist,  and  a normal 
to  the  middle  surface  shall  remain  normal  and  unstrained  dur- 
ing deformation.  These  assumptions  restrict  the  applicability 
of  the  system  of  equations  to  thin  shells  (h/a  <_  20)  . 

A. 2 Derivation  of  Equations  of  Motion 

A polar  coordinate  system  will  be  used  in  the  deriva- 
tion of  the  equations  of  motion.  The  typical  shell  element  will 
be  formed  by  cutting  the  shell  with  two  meridians  having  an  angle 
d8  between  them  and  two  parallel  circles,  as  shown  in  Figs.  1 
and  2. 


s 
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The  six  equilibrium  equations  will  be  determined 
by  summing  forces  and  moments  in  the  x,  y,  and  z directions 
as  shown  in  Fig.  2.  The  z direction  is  normal  to  the  surface 
of  the  element  at  the  center  of  the  element.  The  y-direction 
is  tangent  to  the  middle  surface  and  parallel  to  the  two 
circles.  The  x-direction  is  tangent  to  the  meridian  that 
equally  divides  the  element  and  is  perpendicular  to  the 
circle. 

Summation  of  forces  in  the  x,  y,  and  z directions 
yields  respectively  in  the  following  force  equilibrium  equa- 
tions : 
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Likewise,  summation  of  moments  in  the  x,  y,  and  z directions 
yields  respectively  the  following  moment  equilibrium  equations 
of  motion: 


. , 3M*0  3M0  2 
asin$  — jg—  + a -gg-  + 2acos$M.Q  - a sin$Q 


90 
h/2 


e 


r 2 2 

- psin<j»  / r z dz  = 0 

_wo  3t 


-h/2 

9M.  9M , a 

a3in*  ♦ a -jf  * 


2 . 


acosif)  (M,  - M Q)  - a sin$Q, 
cp  □ ( p 


h/2 


- psinsj) 


/« 


2z  dz  = 0 

at 


-h/2 

h/2  2 

- psinsf)  /r2z  — ^ dz  = 0 


in$  yi 


-h/2 


at 


(A. 2) 


It  has  been  assumed  that  no  external  moments  have  been  applied 
to  the  shell. 


A. 3 Determination  of  Stress  Resultants  in  Terms  of  Displacement 
The  stress  resultants  of  a spherical  shell  as  shown  in 
Figs.  1 and  2,  are  defined  in  terms  of  the  stresses  of  a shell 
segment  as  follows: 


i i 
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- h/2 


- * / °9  ( dz 


“ hi/ 2 


N*e  " / %6  ( d*  - "e 


-h/2 


*♦  " ' / % ( ¥)  2 d2 
t Vi  *\ 


~h/2 


H*e  - m6 * * * / T»e  ( Sr)  2 d2 

-h/2 


/ v 

-h/2  N 

h/2 

=6  - / Te  (*“) «» 

-h/2 


(A. 3) 


In  order  to  relate  the  stress  resultants  and  the  displacement 
quantities,  the  elastic  stress-strain  law  is  used  to  obtain 
stresses  in  terms  of  strains,  which  are  themselves  written 
in  terms  of  displacements. 
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The  strains  can  be  derived  by  considering  an  element 
of  the  shell  as  in  Figure  5 which  describes  the  deformation 

associated  with  the  meridian  strain  e,.  The  undeformed 

<P 

length  of  the  arc  AC  is  ds . Upon  deformation  the 
length  of  the  arc  is  increased  proportionately  with  the  in- 
crease in  the  length  of  the  radius 

(ds  + 30  d$)  T" 


The  change  in  the  length  of  the  arc  A ds  is 

(ds  + d<|>)  1 + J -ds 


or 


. w . 3u  , 
ds  x * 3?  d* 

Therefore , 


_ Ads  w 3u  1 

c<f>  ds  r 3p  r 


Similarly,  it  is  found  that 


IMHHHHMHMWMMBMM 
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u 


= — cot  <J>  + 


3v  w 


e 


r 


3w 

Ir 


1 3u  1 

r sin  <J>  30  r 


V cot  <t> ) 


Y 


0r 


3v  _ v + 1 3w 

3r  r r sin  <f>  30 


3u  _ u , 1 3w 
3r  r r 3 ^ 


Since  we  are  dealing  only  with  thin  shells  and  have  stated  that 
Love's  first  approximation  applies,  a condition  of  strain 
exists  in  which 


0r 


0 


and  w does  not  change  with  r.  Hence,  the  expression  for 
y0r  and  y reduce  to  the  following: 

Y _ 9v  v + 1 3w  _ 0 

0r  ~ 3r  ” r r sin  <J> 


Y ^ 3u_  u ( 1 3w 
$r  3r  ” r r 34> 
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The  above  two  partial  differential  equations  can  be  solved 
for  u and  v with  the  condition  that  u=u  and  v=v  at 
r=a.  The  solutions  are 


u 


a+z 

a 


v 


W.  z 


a sxn  <Ji 


v 


a+z 

a 


Substituting  these  values  back  into  the  three  non-zero  strain 
equations  yields  the  following: 


= ” + + u «+*] 
r 3 <t>  |_  a a ! 

a?z  + aTz  [ ' + % I4!1)] 

“ 5V?[W  - + ^ <a+2>] 


Similarly, 


.-Lf 

a+z  [ 


w - cot  4> 
<J>  a 


weo 2 


a sin  <j> 


t~  + rrin  ♦ (a+z) 


r ■ l ■ r 


I 
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The  shell  is  now  assumed  to  be  in  a condition  of 


I 

I 


2z 

a+z 


w. 


a sin  <J> 


Wq  cos  <p 

— tt7 

a sin  <p 


ua  va  v 

® r + — i - — cot  <j) 

a sm  <J>  a a 

(A. 4) 
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The  remaining  stresses  of  the  shell  can  be  expressed  in  terms 
of  the  strains  as  follows: 


’♦  “ ^7  (e<t>  + VEe> 


°e  = " ? (e6  + vV 


l-v 


T$0  2 (1+v)  Y$0 


(A. 5) 


Substitution  of  Eqs.  (A.  4)  and  (A.  5)  into  Eqs.  (A.  3)  yields 
the  stress  resultants  in  terms  of  the  displacements  of 

the  middle  surface  of  the  shell. 


Na  = 


Eh  1 

w + u.  + v 

a (l-v  ) «■  ♦ 


(w  + sliT^  + u cot  .♦>; 


Ne  - 


Eh  r , v9 
*—  W + —I z 

a(l-v2)  sin  ♦ 


+ cot  4>  u + v(w  + 


v] 


N 


Eh 


r 

[_  a in 


<p6  2a  (1+v)  sin  <p  <p 


+ V . - v cot  <J> 


Eh 

12a2  (1 


) r 

7~ 

.-v2)  ' 


- U . + V (w . cot  <b  + 
<p 


60 


sin 


2"  <(>  ” sin  <J> 


u cot 
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APPENDIX  B 


REDUCTION  OF  ORDER  OF  THE  SHELL 
EQUATIONS  OF  MOTION 


The  reduction  of  the  order  of  the  shell  equations 
is  a necessary  step  in  the  process  of  numerically  determining 
independent  solutions  to  the  equations.  The  eighth  order 
shell  equations  of  motion  can  be  reduced  to  seventh  order 
if  one  solution  is  known.  In  this  case,  the  known  solution 
is  obtained  by  straightforward  numerical  integration  of  the 
differential  equation. 

Let 

u2  = (Uj  + c>  J y2dx 

v2  * V1  /yl  + *2 
W2  * ul/yl  + *3 
s2  - sl/yl  + y4 
N«2  " "«/  yl  + y5 

M*2  ” "u/  yl  + y6 


> 


I 

I 

I 

I 

I 


[[ 

n 

n 
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(sin  <fi  V^)2  = (sin  <t>  Yjdx  + Y-j 


(sin  <(>  V0)2  = (sin  $ VQ)  ^ ^ 


dx  + y 


8 


Differentiation  of  the  above  equations  yields: 


u2  = (u 


2 + c)yL  + ux J y2dx 

v2  - viyi  * vi  Jy id* + y'2 
w2 ■ uiyi + wi  / + yj 

iyi  + a’Jy idx  * 


s2  - 8 


(B.l) 


"♦2  ■ N*  yi ♦ \ J yf*  * y's 

"h  - *Vi + % Jy idx  + y'6 

(sin  0 V^)  2 » (sin  (j>  V^)  ^y^  + (sin  «)>  V^)  [Jy^x  + Y-, 
(sin  <$>  VQ ) 2 = (sin  <J>  Vg)^  + (sin  <J>  VQ)  iJy1 


dx  + yn  (B . 2) 
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Substituting  Eqs.  (B.l)  and  (B.2)  into  Eqs.  (4.5)  and  simpli- 
fying, we  get  the  following  system  of  equations: 


(ux  + c ) y x + -C°— --  c 


hi- 


vn 


a sm  p *2 


Y?  ~ 


(1+v) 


y3  + iF  = 0 


Yo  + vy0  - 
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2 a sin 


:/yi dx 


cot  4> 


y2  " 7? 


2nD 


a E sin  p |_ 


3 cot  p 


- y, 


* a _ n 

E ' k (1-v)  sin  p 


y3  - y4  + wlYl  = 0 


y4  + Siyl  - 


vn 


(1+v) 


2 . 2 
a sin 


y3  + 


y c yfi 

+ + -i  = o 

aE'  D 


y 5 + N*  yl  - "7 
1 a 


(1-v2)  cot2  <{> 


cfy .dx  - <i-v2)y 

1 sin  <j> 


(l-v2)E'  cot  <j>  2n2D  cot  P ,,  ..v 

2 ~T — 2~~  (1_v) 

a k sin  p 


y3  ~ .3 


2n  D(l-v) 


a k sin  <j> 


+ HryJ  JEot  <P  v + ^ 


-hi* 


s + s-inr?  — rr:  (k'2)y 

r a sin  p 


- ( phui  )c  / y,dx.  * 0 

.2, 


y’  ♦ H^Yl  + -(--XV0t  ♦ cj Vjdx  OOt  ♦ + JXFT  y2 
+ nW^c^(|+1  + v y3).Uq»D 

a Sill  <P  3t 


2n 


(_k  sin  p 


+ (1+v)  cot  <)> 


y4  ♦ «i-v,  ye  * nrj 


n 

a 
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2nDy8 

~2  • 7~T 

a E k sxn 


= 0 


y;  + (sin  ♦ v^lyj  + ll~sIE 


,,  . n D ( 1+v) 

(1  V)  + ~2P"  _,_2 


a E sin  <j>  J 


:Pi‘ 


(c  / yndx  cos  4>  + ny_)  + ■ ^ E- 

a sin  4> 


. 2 

sin  <p  + 


4 2 2 

n D . 2n  D cot  4> 

— 5 — i — ' *y~  1 + — v'  , ' ~ r' 

a e sinz  <()  azEk(l+v) 


(l-v)E  cot  4> 

1 Q a 


n (l-v)D  _ 2nD 

2 I ” “2  ■ 

a E sin  <j>  a E k sin  <)> 


(1+v)  sin  » 
y4  a . 


+ Vs Y*  + 2,nD  Cot  ♦ y.  + Pho.2  sin  ♦ y,  - 0 

a sin  <(>  aJE'k  sin  <p 


Yg  + (sin  <p  VQ)y1  - 


nE 


cot  <t>  (l-v2)k  c J y2dx  + -gg  ^(l-v2)ky2 


- Jr  [E,(l-v2)(l 


H*v-  y + V 
IT  y5  r y6 


l n2° 

“2 — ~T~ 
a sin  <(> 


)I 


D ( 1+v)  cot  <t> 
naE ' ^ 


" iPsfnS  y8  + (phw2  Sin  *)y2  " 0 

(B.3) 


Hence,  the  order  of  the  original  system  of  shell  equations  has 
been  reduced  from  eight  to  seven.  This  seventh-order  system 
can  now  be  numerically  integrated  to  obtain  non-growing  solu- 
tions as  described  in  section  2.3. 
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APPENDIX  C 


SHELL  DIFFERENTIAL  EQUATIONS  OF  MOTION  WRITTEN 
IN  TERMS  OF  THE  INDEPENDENT  VARIABLE  Z 


In  the  shell  equations  the  independent  variable 
<t>  is  the  colatitude,  which  ranges  from  0°  to  90°.  Dif- 
ficulties are  encountered  when  trying  to  integrate  the  shell 
equations  numerically  in  the  neighborhood  of  the  apex.  This 
is  due  to  the  fact  that  many  of  the  coefficients  in  the  equa- 
tions are  of  the  form  c/sin  <(>,  which  gets  very  large  as 
4>  approaches  zero  and  is  undefined  at  <(>=0.  If  the  indepen- 
dent variable  in  these  equations  is  changed  from  <$>  to 
z=ln(sin  <f>)  this  problem  can  be  eliminated. 

If 

z = ln(sin  4>)  (C.l) 


then 


dz  m cos  <t> 

d<t>  sin  <j)  (C.2) 


Consider  Eqs.  (4.5a) 


1 dw 
a cT$ 


s 


- 


! 


can  be  expressed  as  follows: 

dw  _ dw  dz 

d<j>  dz  ' d4>  (C.  3) 

Substitution  of  Eqs.  (C.2)  and  (C.3)  into  Eqs.  (4.5a) 
yields 


_1  dw  cos  4) 
a dz  * sin  <p 


or 


1 dw  sin  <p 
a dz  cos  $ 


Similarly,  the  other  equations  in  terms  of  z are  written  as 
follows : 


2 

ds  = / vn 

adz  = v^a  cos  <J>  sin  4> 


(i*vyn  f) 


du  _ _ v _ vn 
ad 0 a u a cos  <J> 


(1+v)  tan  <t> 


v - 


a 
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dv  _ n + v + 2nD 

ad<(.  a cos  4>  a a2£l  CQS\ 


cot  <{> 

a 


' s) 


. 2 tan  4>  „ 

+ E'kd'-v]  V8 


= !i  [(!-*»>  cot  ♦]  U 


+ -5-—' (1~V2)V  + 


a sin  4> 


(1-v2)E’ 


. 2ff2q(l-W  „ + 2n2D(l-v) s - ilTiL  tan  * 
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« - U--^)P  (u  cot  * t " N . n2D(l-W  (2  + , + v)  „ 

ad<p  p"  V sm  4>y  a3  sin2  ^k  ; 
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a cos  4> 


4 2 2 
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— * 5 — + — * w 

a*E'  sin  4)  a E’k(l+v) 


(l-v ) E * n* ( 1-v ) D 


2n2D 


a E'  sin  0 a E'k  sin 


s (1+v)  sin  <»  „ 
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~2 — M,  + -2^D — Vg  + pho)2  tan  4>  sin  4>  w 

a*  cos  4>  ^ aJE'k 


53.  ve  sin  ♦ m^[  (1_v2)ku  + (1-y2)kv  + TT4 

e.(1_d2,  ( 1 t ^"2d;  )]w  + s 

'a  sin  4>/  J 
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na  v tan  <J>  „ v tan  <t>  .. 

£n “ — ttt — *■  M 

E <P  E < 


_ sin  <p 

a v6 


(phw  tan  4»  sin  <j>)v 


(C.4) 
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APPENDIX  D 

DETERMINATION  OF  THE  INITIAL  EIGENVECTOR 
IN  THE  ROBINSON-HARRIS  METHOD 

The  initial  eigenvector  used  in  the  Robinson-Har- 
ris  method  is  determined  from  the  latest  initial-value  sol- 
utions in  the  Holzer  method.  These  initial  value  solutions 
are  linearly  combined  to  form  a solution  which  will  come 
close  to,  but  not  completely,  satisfy  the  homogeneous  boun- 
dary conditions.  This  solution,  in  turn,  becomes  the  first 
approximation  of  the  eigenfunction  in  the  Robinson-Harris 
method. 

The  problem  to  be  solved  in  determining  the  com- 
plete solution  is  of  the  form 


DX  = T (D.l) 

where  D is  the  matrix  of  the  values  of  the  displacements 
and  stress  resultants,  X is  the  vector  of  unknown  weighting 
coefficients  and  T is  the  complete  solution.  To  deter- 
mine X we  evaluate  Equation  (D.l)  at  the  boundaries,  that 
is. 


DX  = 0 


(D.  2) 
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or  for  n=0 


u1(0) 

u2(0) 

u3(0) 

c 

itk 

0 

c 

u> 

0 

u6(0) 

xi 

s^O) 

s2(0) 

s3  (0) 

s4  (0) 

s5(0) 

s6(°) 

x2 

o 

-e- 

H * 

O 

Q*  <°> 

Q*  (°) 
♦3 

Q* 

V01 

Q.  <°> 
*6 

x 

x3 

w 

u2{<*)o) 

U3(V 

U4(V 

w 

W 

X4 

w 

W2(V 

W3(V 

W 

*5 ‘♦o’ 

w 

X5 

S1(V 

S2(V 

s3  ^o) 

w 

*5(V 

W 

x6 

(similar  equations  can  be  written  for  n=l  and  n>l) 

(D.3) 


If  we  assume  that  x^=-l  then  we  can  formulate  a system  of  five 
equations  in  five  unknowns  in  terms  of  the  remaining  x's. 
Performing  the  indicated  multiplication  in  Eqs . (D.3) 
yields 


u3  (OJXj^ 

4*  • • 

• • 4* 

V0)X6 

81(0)xl 

4“  • • 

• • 4- 

S6(0)X6 

Q.  (°)Xi 

4-  • • 

• • 4- 

ul(*o,xl 

+ • . 

• • 4* 

ul(*o)x 

W1(^0)X1 

+ • • 

• • 4- 

wl(*o,xi 

S1(*0)X1 

4*  . . 

• • 4- 

sl^o,xi 
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Rearranging  terms  we  get 


u2(0) 

. . . . Ug(0) 

ux(0) 

O 

CM 

(0 

. . . . s6(0) 

X2 

sl(0)xl 

s1(0) 

Qa  (0) 

....  Q.  (0) 

X3 
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%,0,xl 

Q*  (0) 
♦l 

u2(*o} 

. . . . u6(<J>o) 

X4 

ul<*o,xX 

u^O) 

W2(V 

* • * * w6t*o) 

x5 

WX1 

w.  (0) 

-L 

S2(V 

• * * * s6(<,*o) 

X6 

=1<VX1 

s1(0) 

(D . 4 ) 


Equation  (D.4)  can  be  written  symbolically  as 


°6x5  x E5x1  = C6xl 


(D.5) 


Pre-multiply  both  sides  of  Equation  (D.5)  by  D transpose 
to  get 

D5x6  x D6x5  x E5x1  " D5xl  X C6xl 


Equation  (D.6)  consists  of  five  equations  in  five  unknowns 
which  can  be  solved  for  Xj,.  • , . ,Xg.  The  weighing  co- 
efficients are  then  substituted  into  Equation  (D.l)  to 
give  us  the  complete  solution,  which  is, 

6 

ih  ui*i  " UT 


e 

i-1  sixi  = ST 


ill  wixi  - "T 


6 

iSi 


N,  xi 

♦i  x 


6 

i-1 


x . 

X 


6 

i£i  ■ °*T 


In  general,  when  the  initial  eigenvector  is  sub- 
stituted back  into  Eqs.  (D.2),  the  boundary  conditions 
will  not  be  satisfied  exactly.  The  residuals  from  these 
equations  will  be  used  to  compute  the  new  distribution 


factors. 
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APPENDIX  E 


DERIVATION  OF  THE  KIRCHHOFF  FORCE 

Kirchhoff's  shear  forces,  represented  by  and 
Vg , will  be  derived  in  this  appendix.  These  forces  are 
caused  by  the  effects  of  the  stress  resultants  N^0,  M^0, 

and  at  the  edge  of  the  shell. 

Figure  6 illustrates  two  adjacent  incremental  seg- 
ments of  the  edge  of  a shell.  The  location  of  the  angles 
a , <p  and  0 are  given  in  Figure  7.  The  moments  in  Figure  6 
can  be  replaced  by  force  equivalents  as  shown  in  Figure  8. 
From  Figure  8 it  can  be  seen,  that  because  of  the  slight 
divergence  there  is  a horizontal  resultant  force  Ffi  sin  d 
acting  on  the  section.  Also  from  Figures  8 and  11  we  see 
that  because  of  the  curvature  of  the  element  ds  we  have  the 
resultant  force  F =sin  4 F d6  (E.2).  Therefore,  the  total 
force  per  unit  length  in  the  direction  of  N^0  (see  Figure  1) 
is  v0=N^e-Ft/ds  (E.3).  Substitution  of  Eqs.  (E.l)  and  (E.2) 
into  (E.3)  yields 


V 


6 


(sin  <|>  F d0)/(a  sin  $ d0) 
n 


V 


0 


- V« 


' 


which  reduces  to 


L35 


From  Figure  C we  see  that  F =M._  Thus, 

n <po 

ve  “ VM*e/a 

A similar  reasoning  may  be  used  to  derive  the 
force  which  acts  on  the  edge  of  a shell  in  the  same 

direction  as  (see  Figure  1) . We  see  from  Figure  11 

that  the  total  shear  force  on  that  section  is 

v»  = ° ♦(frd“)/ds 

or 

V ° +(tt  d0)'ds  (E 

Substituting  Eqs.  (E.l)  into  (E.4)  we  get 

V - Q$  ♦ TT  (d8/a  sin  ♦ 

From  Figure  8 we  see  that  F =M.„.  Therefore, 

n 

„ _ . . 3M<t>e 

V<)>  + a sin  $3  0 

Using  Eqs.  (4.1)  from  Section  4.2  we  get 
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Differential  equations;  numerical  integration;  edge  effects 
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An  accurate  procedure  Is  described  for  numerically  solving  two-point 
boundary  value  problems  which  contain  growlnq  solutions.  The  procedure 
Involves  the  process  of  reducing  the  order  of  a differential  equation  when 
one  solution  Is  known.  Two  applications  of  the  procedure  are  given,  a 
fourth  order  differential  equation  with  two  growing  solutions  and  a system 
of  eighth  order  differential  equations  of  motion  for  a hemispherical  shell. 
In  both  examples  before  the  procedure  Is  started,  the  equations  are  re-  — ; 
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written  as  a system  of  first  order  differential  equations.  It  was  found 
that  when  solving  two-point  boundary  value  problems  by  the  reduction  of  order 
method,  first  order  differential  equations  were  generally  easier  to  work 
with  than  higher  order  differential  equations.  For  both  applications  a 
computer  program  was  developed  to  solve  the  system  of  differential  equa- 
tions. 
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